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Abstract. The multiple scattering of scalar waves in diffusive media is investigated by 
means of the radiative transfer equation. This approach amounts to a resummation of the 
ladder diagrams of the Born series; it does not rely on the diffusion approximation. Quan- 
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within the diffusion approximation. 
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1. INTRODUCTION 



The theory of multiple light scattering has been a classical subject of interest for 
one century, which attracted the attention of many scientists, including Lord Rayleigh, 
Schwarzschild, and Chandrasekhar. Standard books are available, such as those by Chan- 
drasekhar [1], Ishimaru [2], and van de Hulst [3]. The discovery of weak-localization effects, 
such as the celebrated enhanced backscattering cone [4] , yielded a revival of theoretical and 
experimental work in the area of multiple scattering in disordered media. Much progress 
has been done recently in the analysis of speckle fluctuations [5], in analogy with the 
conductance fluctuations observed in mesoscopic electronic systems [6]. 

Laboratory experiments are often performed either on solid Ti02 (white paint) sam- 
ples, or on suspensions of polystyrene spheres or of Ti02 grains in fluids. In most cases, 
the wavelength Ao of light in the diffusive medium, the scattering mean free path £, and 
the thickness L of the sample obey the inequalities Ao <C £ <C L. Multiple scattering 
is also of interest in biophysics and medical physics, in order to understand the trans- 
port of radiation through human and animal tissues. Besides light, all kinds of classical 
waves undergo multiple scattering in media with a high enough level of disorder, i.e., of 
inhomogeneity. Well-known examples are acoustic and seismic waves. The propagation 
of electrons in disordered solids also pertains to this area, since Quantum Mechanics also 
basically consists in wave propagation. Multiple scattering of waves in disordered media 
admits the following three levels of theoretical description: 

(i) The macroscopic approach consists in an effective diffusion equation, which describes 
the transport of the diffuse (incoherent) intensity /(r, t) at point r at time t. This approx- 
imation turns out to be very accurate in the bulk of a turbid medium, and more generally 
on length scales much larger than the mean free path I. The diffusion approximation 
yields several interesting predictions, among which we mention the 1/L-decay of the total 
transmission through an optically thick slab of thickness L » f, the decay time of the 
transient response to an incident light pulse, or the memory of a typical speckle pattern 
when the frequency of light is varied. This approximation also allows for a quantitative 
prediction of the diffuse image of a small object in transmission [7] . 

(ii) The mesoscopic approach, used by astrophysicists throughout the classical era of the 
subject, is referred to as radiative transfer theory [1-3]. This theory relies on the radiative 
transfer equation, which is a local balance equation, similar to the Boltzmann equation in 
kinetic theory, for the diffuse intensity i"(r, n, t), with n being the direction of propagation. 
This approach leads in a natural way to distinguish between the scattering mean free path 
£ and the transport mean free path £* , to be defined below. The diffusion approach (i) is 
recovered in the limit of length scales much larger than £* . 

(iii) The microscopic approach consists in expanding the solution of the wave equation in 
the disordered medium in the form of a diagrammatic Born series. In the regime £ ^> Ao 
the leading diagrams can be identified, in analogy with e.g. the theory of disordered su- 
perconductors [8] . For the diffuse intensity they are the celebrated ladder diagrams, which 
are built up by pairing one retarded and one advanced propagators following the same 
path through the disordered sample, i.e., the same ordered sequence of scattering events. 
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This picture agrees with that behind the radiative transfer equation, which is a classical 
transport equation for the intensity. The ladder diagrams can be resummed and yield an 
integral equation of the Bethe-Salpeter type for the diffuse intensity, which we refer to as 
the Schwarzschild-Milne equation. The radiative transfer approach (ii) is thus recovered 
in the weak-disorder (£ 3> Ao) regime. A further step consists in including some of the 
subleading diagrams, of relative order Xq/£ <C 1, which account for interference effects be- 
tween diffusive paths. Among those contributions, the class of maximally-crossed diagrams 
is of particular interest, since it describes the aforementioned enhanced backscattering phe- 
nomenon [4]. 

To summarize this discussion, each of the descriptions mentioned above represents a 
qualitative improvement with respect to the previous ones (see e.g. ref. [9]). In order to 
derive quantitative estimates of physical observables in the regime i Ao of experimental 
interest, it is sufficient to consider radiative transfer theory. The macroscopic approach (i) 
is clearly insufficient, since we aim, among other things, at a description of the crossover 
from free propagation to diffusive transport which takes place in the skin layers of thickness 
of order £, near the boundaries of the disordered medium. This phenomenon is, by its very 
definition, beyond the scope of the diffusion approach. The radiative transfer approach 
(ii) leads to study the Schwarzschild-Milne equation (2.24). This equation has only been 
solved exactly in a limited number of cases. Analytical results have been obtained for 
isotropic scattering of scalar waves [1, 10-12], and in the case where the scattering cross- 
section depends linearly on the cosine of the scattering angle [1, 13]. The case of anisotropic 
scattering has essentially been investigated by means of general formalism and by numerical 
methods [2, 3, 14]. 

For a finite sample, physical observables such as the reflected or transmitted intensity 
in a given direction depend on the particular realization of the sample under consideration. 
Such quantities are indeed the results of intricate interference patterns throughout the 
sample; they only become self-averaging quantities in the limit of a large enough sample. 
This definition can be made precise by means of the dimensionless conductance g ~ Af£/L, 
related to the number Af ~ A/X^ of open channels, where A is the transverse area of the 
sample. The self-averaging regime corresponds to g 3> 1. The whole distribution of 
observables is therefore of physical interest, as long as g is not very large. We mention 
a recent experiment [15], where the third cumulant of the total transmission, an effect 
of relative order l/g 2 , has been measured and compared to theoretical predictions. This 
third cumulant is of the same order of magnitude as the universal conductance fluctuations, 
either in electronic [6] or in optical [5] systems. In fact the full distribution of the total 
transmission through an optically thick slab has been derived recently [16]. In the following 
we focus our attention on the mean values of physical quantities. 

The vector character of electromagnetic waves also introduces its own intricacy. In 
general four coupled integral equations for two transfer functions have to be solved, which 
are associated with the Stokes parameters of the diffuse light in the medium. These equa- 
tions have been solved exactly in the case of Rayleigh scattering, i.e., the regime where the 
size of the scatterers is much smaller than the wavelength [1] . Among specific features per- 
taining to diffusive light propagation, let us mention the dependence of the backscattering 
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enhancement factor on the polarization state of the outgoing diffuse radiation, for a linearly 
polarized incident beam [17, 18], or the progressive destruction of the backscattering peak 
induced by a magnetic field, due to the Faraday rotation in a magneto-optically active ma- 
terial [19, 20]. The latter effect has been recently observed and interpreted quantitatively 
[21]. ' 

Furthermore, in practical situations the optical index no of the scattering medium is 
often different from the index n\ of the surrounding medium. This index mismatch, mea- 
sured by the ratio m = n$/ni, causes reflections at the interfaces. In the regime of a large 
index mismatch (m <C 1 or m > 1), the transmission across the interfaces is very small, 
so that the light is reinjected many times in the diffusive medium. As a consequence, the 
skin layers become very thick in this regime. More generally, the diffusion approximation 
works better and better as the index mismatch gets large. Improvements of the diffusion 
equation have been proposed [22, 23], which take internal reflections into account. It is in 
fact possible to derive analytical expressions for the reflected and transmitted intensities 
and other observables in this regime. This asymptotic analysis has been performed in 
ref. [12] for isotropic scattering. It will be generalized hereafter to the case of general 
anisotropic scattering. This approach provides accurate results, even for a moderate index 
mismatch m. 

More generally, one of the main goals of this paper is to quantify the dependence 
of physical quantities on the anisotropy of the scattering mechanism. It has been known 
for long from radiative transfer theory that two length scales are involved in the case 
of general anisotropic scattering: the scattering mean free path £, namely the effective 
distance between two successive scattering events, and the transport mean free path £* , 
namely the distance over which radiation looses memory of its direction. Both mean 
free paths, to be defined more precisely in section 2, depend on details of the scattering 
mechanism, such as the shape, size, and dielectric constant of the scatterers. In the regime 
of very anisotropic scattering, we have £* ^> £. The ratio r* = £* /£ ^> 1 will be referred 
to as the anisotropy parameter. In some experimental situations r* can be of order ten or 
larger [24]. The regime of most physical interest is then A < f < f < L. 

The setup of this paper is as follows. In section 2 we present some general formalism on 
the radiative transfer equation, and we derive the associated Schwarzschild-Milne equation. 
We show how solutions of the latter equation yield predictions for quantities of interest, 
such as the diffuse reflected and transmitted angle-resolved intensities, and the shape of 
the enhanced backscattering cone. Section 3 is devoted to the large index mismatch regime 
for general anisotropy. In section 4 we derive a complete analytical solution of the problem 
in the regime of very anisotropic scattering (i.e., a cross-section strongly peaked in the 
forward direction), in the absence of internal reflections. The paper closes up with a 
discussion in section 5. 

2. GENERALITIES ON ANISOTROPIC MULTIPLE SCATTERING 

Throughout the following we restrict the analysis to multiple scattering of scalar waves 
by scatterers located at uncorrelated random positions, in the regime where the scattering 
mean free path £ is much larger than the wavelength Aq of radiation in the medium. We 
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have summarized some useful notations and definitions in Table 1. 

As stated in section 1, we put special emphasis on the dependence of physical quanti- 
ties on the anisotropy of the scattering mechanism. It is known that, after averaging over 
the random orientations of the individual scatterers, the differential scattering cross-section 
of arbitrary anisotropic scatterers can be written in the following form [1, 2] 

da(n, n') = (u/4n) 2 p(G)dn' . (2.1) 

In this formula u is the scattering length, n and n' are unit vectors in the incident and 
outgoing directions, 6 is the angle between these directions, so that cos© = n.n', and dfl' 
is an element of solid angle around the direction n'. 

We assume that there is neither absorption, nor inelastic scattering, i.e., that the 
albedo is unity [the situation of non-conservative scattering will only be considered in 
section 2.6]. The phase function p(Q) then obeys the normalization condition 

f dO' , . f 1 dcosG 1 

J 4^(e) = — P(e) = i. (2.2) 

The total cross-section reads a = u 2 /(4tt), and the scattering mean free path £ is given by 

1=- = ^ (2.3) 
no nu z 

where the density n of scatterers is assumed to be small, in such a way that we have 
A . 

In the particular case of isotropic scattering, the phase function p(6) = 1 is constant. 
In the general case of anisotropic scattering, the phase function p(0) is a non-trivial 
function, often peaked in the forward direction, albeit in a more or less prominent way. 
As recalled in section 1, one has to distinguish between the scattering mean free path £, 
which is the typical distance between two successive scattering events, and the transport 
mean free path £*, which represents the distance over which radiation looses memory of its 
direction. The dimensionless ratio of both mean free paths, 

£* 1 

r * = J = i - (cos e> ' (2 ' 4) 

will be referred to as the anisotropy parameter. We have introduced the notation 

/ ^^ f 1 dcosG _ , n 

(cos0)= / cos0p(e). (2.5) 

J-i 2 

We commonly have £* > £, since scattering usually preferentially takes place in the forward 
direction. The regime of very anisotropic scattering, where p(@) is strongly peaked around 
the forward direction, corresponds to £* 3> £. This regime will be investigated in detail in 
section 4. 
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2.1. General formalism and Schwarzschild-Milne equation 



In this section we present some general formalism on anisotropic multiple scattering, 
extending thus the treatment of the isotropic case presented in ref. [12]. Some of the 
results exposed below are already present in lecture notes by one of us [9]. 

We begin with a derivation of the Schwarzschild-Milne equation. This key equation 
of radiative transfer theory will be the starting point of our analysis. In the regime I ^> Xq 
under consideration, in the absence of internal sources of radiation, and in stationary con- 
ditions, the quantity of interest is the specific intensity J(r, n) of radiation at the position r, 
propagating in the direction n. The specific intensity obeys the time-independent radiative 
transfer equation, that takes the following local form 



is commonly referred to as the source function. We recall that the phase function p(n, n') = 
p(0) only depends on the scattering angle 6. 

As recalled in section 1, the radiative transfer equation (2.6) [1-3] can be considered as 
a mesoscopic balance equation for the light intensity inside the diffusive medium, somewhat 
analogous to the Boltzmann equation in the kinetic theory of gases. It is equivalent to 
the Bethe-Salpeter equation, obtained by resumming the ladder diagrams of the Born 
series expansion of the intensity Green's function of the problem. These diagrams are the 
dominant ones for £ 3> Ao, i.e., to leading order in the density n of scatterers. 

We consider a sample of diffusive medium in the form of a slab of thickness L, limited 
by the two parallel planes z = and z = L. The mean optical index no of the slab can be 
different from the index n\ of the surrounding medium. The index mismatch, measured 
by the ratio m = no/ni, generates internal reflections at the interfaces. We introduce the 
optical depth r = z/t of a point in the sample, and the optical thickness b = L/£ of the 
sample. Most observables of interest can be derived by integrating out the (^-dependence of 
the intensity 7(r, n), i.e., by only considering quantities with cylindrical symmetry around 
the normal to the slab. The radiative transfer equation (2.6) then takes the form 



£n.V/(r, n) = T(r, n) - J(r, n). 



(2.6) 



The quantity 




(2.7) 



(2.8) 



or equivalent ly 



^[/(r^)e^] =r(T,p)e T /". 



(2.9) 



The source function T(r, y) is related to the specific intensity 7(r, y) through 



T(T,y) = V[l(r, fJ ,)], 



(2.10) 
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where we have introduced the integral operator V 

•1 A.. I 



*W/0] = J J^-PofaSW), (2-11) 
which involves the reduced phase function 

Po(/j,, fjf) = J ~^ p ^ <P, V', V 7 ')- (2-12) 

The original phase function p(G) and the reduced one po(/^, fj/) are more simply related 
[1] through their expansions in the Legendre polynomials 

p(G) = J](2£+l)^P £ (cose), 

e ~° (2-13) 

£>0 

We mention for further reference that the Legendre polynomials |P^(/x),£ > 0} form 
a complete set of orthogonal eigenfunctions of the integral operator V, introduced in eq. 
(2.11), with eigenvalues zj£, i.e., 

V[P e (jj)] = w e P £ (fi). (2.14) 

The eigenvalue wq = 1, associated with the constant eigenfunction Pq(h) = 1, repre- 
sents the unit albedo of the problem. More explicitly, 

jf ^Po(/i,/*') = l. (2.15) 

This identity yields the zero mode of the diffusion equation in the long-distance limit. 

The first non-trivial eigenvalue w\ is also of special interest. It is indeed related to 
the anisotropy parameter r* of eq. (2.4). Since P\{p) = £t, we have 

Wl = (cos0) = 1 - (2.16) 

More explicitly, 

d ' 

— po(fi, n')n' = win = (cos Q)/j,. (2.17) 

We also notice that all the eigenvalues of the operator V are trivial in the particular 
case of isotropic scattering, since wi = for I > 1. 



For simplicity, we consider first a half-space geometry (6 = oo). We assume that the 
limiting plane r = of the sample is subjected to a cylindrically symmetric incident beam, 
characterized by an angle of incidence 9 a , i.e., that the incident intensity does not depend 
on the azimuthal angle <p a . This is automatically satisfied at normal incidence (9 a = 0). 
Under these circumstances, the inward intensity on the limiting plane r = + contains 
both the normalized refracted incident beam, coming in a direction defined by fi a , and the 
intensity coming from the bulk of the medium, after being reflected once at the interface. 
Using the radiative transfer equation (2.8, 9), we thus get 

I(0+,fi,fi a ) = 2S(fi-fi a ) + ^- / dre-^rXr,-^) (// > 0), (2.18a) 

A* Jo 

where the Fresnel reflection coefficient R(fi) is given in Table 1. The only contribution to 
the outward intensity on this plane comes from the light rays which have experienced at 
least one scattering event, namely 

1 f°° 

I(0 + ,-H,Ha) = - dre- T ^T(r,-fi,iJ a ) (// > 0). (2.18b) 
A* Jo 



The radiative transfer equation (2.8, 9), together with the boundary conditions (2.18) 
at t = + , can then be recast in the integral form 

I(r, fi, fi a ) = 28{n - fi a )e~ T /^ + (K*T) (r, fi, fi a ). (2.19) 

Here and throughout the following, the star denotes the convolution product 

(K * T)(r, fi, n a ) = J™ dr' J ^ ^-K(t, fi, r\ ^)T(r' fi', fji a ). (2.20) 

The kernel K involved in this integral can be split into two components: K = Kb + K^. 
The bulk kernel Kb contains the contributions to the intensity at depth r arising from the 
scattering from either smaller or larger depths. The layer kernel Kl takes into account 
the intensity being scattered at depth r' in the direction of the wall, then reflected there, 
and then scattered at depth r. These kernels read explicitly 

K B (r, h, t', h') = 25(v - h')6{h{t - T '))±-e-^')/^ 

(2-21) 

K L (r, ^ r', „') = 25(p + ^)^> e -^')/^ 



The final step consists in using eq. (2.10) in order to derive from eq. (2.19) the 
following closed integral equation for the source function 

T(r, fi, aO = po(/x, Va)e- T/ ^ + (M * T)(r, fi, fi a ), (2.22) 
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which we refer to as the Schwarzschild-Milne equation of the problem. 

The kernel M has the following two components, in analogy with the kernel K from 
which it derives 



(2.23) 



(2.24) 



M B (r, r', „') = 9( P '(r - ^fAt^l^-r')/ ^ 

I A* I 

M L (r, ^ r', J) = g (V) ^(^-^) jR(V)e -( T+ ^)/lM'l 

IA 4 I 

so that eq. (2.22) reads explicitly 

r(r, (jl, fJL a ) =p (p,fi a )e- T/ ^ 

+ f dr ' f %P^^~ {T ~ T,)/ "'nr'^'^a) 
Jo Jo Z A* 

+ jT dr' j X ^jpo^ -^)e- (T '- T)/ ^r(r', V, ^) 

+ fdr' f^R{»')po{^»')e-^')/»'T{ T \-Li'^a). 
Jo Jo z f Jl 

In the case of isotropic scattering, the phase function po(fJ>, £*') = 1 is a constant, and 
the source function T(r, /U, jii a ) does not depend on fx. The Schwarzschild-Milne equation 
(2.24) thus takes a simpler form, that has been extensively studied in ref. [12]. The rest 
of this section is devoted to an extension of the results derived there to the general case of 
anisotropic scattering. 

2.2. Solutions of Schwarzschild-Milne equation and sum rules 

In the general case of conservative anisotropic scattering, the Schwarzschild-Milne 
equation (2.24) has a special solution Ts(t, /j,, fi a ) which remains bounded as r — > oo, 
whereas the associated homogeneous equation has a linearly growing solution r#(r, /S). 
More precisely, it can be checked by means of eq. (2.17) that both source functions and 
the associated specific intensities have the following asymptotic behavior for large depths 
(r ^> 1), up to exponentially small corrections 

T S (r, fl, fU a ) ~ Ti(fJL a ), 
I S (r,fi,fi a ) « n(ji a ), 

r H (T,iu)^r-iu(T* -1) + t t*, 
Ih{t, a*) ~ t - /IT* + TqT*. 



The quantities To and Ti(fi a ) are unknown so far. Let us anticipate that they play a 
central role in the following, in the sense that they will bear the full non-trivial dependence 
of physical quantities on the scattering mechanism. These quantities also obey remarkable 
sum rules, to be derived below. 
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To do so, it is most convenient to introduce the Green's function Gs(t, fi, t', fi!) of the 
problem, along the lines of ref. [12]. It is defined as the solution which remains bounded 
as t — > oo of the equation 

G s (r, ijl, r', fi') = po(//, li')8(r - r') + (M * G s )(r, fi, r', fi'). (2.26) 

The kernels K and M and the Green's function Gs possess the symmetry properties 

K(t, h, t', fjf) = K(t', -//, r, -//), 

M(r, ^, r', ^) = M(r', -//, r, -/x), (2.27) 
G s (r, fi, t', fx') = G s (t', -fi', t, -fi), 

which merely express the time-reversal symmetry of any sequence of scattering events. 

As a consequence of eq. (2.22), the special solution r s(r, fi, fi a ) can be expressed in 
terms of the Green's function as 

/•OO 

T s (r, fi, fi a ) = / dr'e- T '/^G s (r, fi, r', fi a ). (2.28) 
Jo 

We also define for further reference the following bistatic coefficient 

"f(fi a ,fib)= dre~ T/lIb Ts(T,-fi b ,fi a ) 
Jo 

rOG POO 

= / dre- T /^ / dT'e- T '/^G s (T,-fi b ,r',fi a ). 
Jo Jo 



(2.29) 



The latter expression defines the bistatic coefficient for any complex values of its arguments 
with Re fi a > 0, Re fib > 0, even outside the physical range fi a , fib < 1- The symmetry 
7 '(a* a, l^b) = lifibiHa) is a consequence of the properties (2.27). It is thus again due to 
time-reversal symmetry. 

On the other hand, a relationship between both solutions Th(t,ijl) and V s(r, fi, fi a ) 
of the Schwarzschild-Milne equation can be derived as follows. The Green's function 
Gs(t, fi, t', fi') is clearly asymptotically proportional to the homogeneous solution Th(t, fi) 
when t' goes to infinity, namely 

Jim G s (r, fi, r', fi') = ^T h (t, fi), (2.30) 

where the proportionality constant D will be shown in a while to be equal to the dimen- 
sionless diffusion coefficient (2.35). 

As a consequence of eqs. (2.25, 28-30), we have 

rioo= is™ ^a,fi b ) = i r dTe -^ rff(T] _ Ma) . (2 . 31) 
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We now turn to the derivation of two groups of sum rules obeyed by the quantities 
defined so far, which are related to the so-called F and K integrals, with the notations of 
Chandrasekhar [1]. 

The first group of two sum rules is a consequence of the conservation of the flux in 
the ^-direction in a non-absorbing medium, given by the following F-integral 

F(r) = J^vJ(T,p). (2.32) 

It can indeed be checked that eqs. (2.8, 15) yield dF/dr = 0. 

We investigate first the F-integral Fs associated with the special solution Ts(r, fj, a ). 
Considering the r — > oo limit determines Fs = 0, whereas considering the r — > limit 
yields the sum rule 

J ^T(p)i{ l i,n a )=n a , (2.33) 

where the transmission coefficient T(/x) = 1 — R(/J,) is given in Table 1. The \i a — > oo limit 
of eq. (2.33), together with eq. (2.31), yields another sum rule, namely 

f^T^)rM = l. (2.34) 

We can evaluate in a similar way the F-integral F H associated with the homogeneous 
solution Th(t, ji). This yields no independent sum rule, but leads to the identification of 
D with the dimensionless diffusion coefficient 

D = r*/3. (2.35) 

The diffusion coefficient indeed reads -D p h ys = c£D = c£*/3 in physical units [1-3]. The 
transport velocity indeed coincides with the velocity of light in vacuum c, to leading order 
in the regime t Ao- 

Besides the sum rules (2.33, 34), which were already given in ref. [12] for isotropic 
scattering, the radiative transfer equation also admits another group of two sum rules, 
which are novel in this context, and whose physical origin is less evident. Consider the 
so-called iif -integral, again with the notations of Chandrasekhar [1] 

K(r) = £J±v?I(t,vl). (2.36) 

It can be checked that eqs. (2.8, 16, 17) yield dK/dr = -F/t*, whence K{r) = -Ft/t* + 
Kq, with K being independent of r. By considering the X-integrals associated with the 
special solution Ts(r, /jl, /j, a ) and with the homogeneous solution Th(t, /jl), we obtain after 
some algebra the following sum rules 

f Q ^(1 + R{n))m^^a) = ^) _ ^ (2.37a) 
I q ^(1 + R(^tM=t . (2.37b) 
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2.3. Diffuse reflected intensity 

The evaluation of the angle-resolved diffuse reflected intensity by means of the general 
formalism exposed above closely follows the lines of ref. [12]. We consider a half-space 
geometry, and we assume that the limiting plane of the sample is subjected to a cylindri- 
cally symmetric incident beam, characterized by an angle of incidence 9 a . This technical 
assumption includes the case of a plane wave at normal incidence (9 a = 0). Under these 
circumstances, the diffuse reflected intensity per solid angle dO& reads 

dR(a^b) R cos 6 a T a T b 

)— '-=A H (e a ,e b ) = - j l{Va,Vb), 2.38 

where we have used some of the definitions of Table 1: m = uq/ui is the index mismatch, 
and T a = T(/z a ) and T b = T(/ib) are the transmission coefficients in the incident and 
outgoing directions, respectively. 

The most non-trivial part of the above result (2.38) consists in the bistatic coefficient 
7(/i a , Hb), whose definition and general properties have been exposed in section 2.2. It will 
be evaluated in the large index mismatch regime in section 3, and in the very anisotropic 
regime in section 4. 

2.4. Diffuse transmitted intensity 

In this section we consider the angle-resolved mean transmission of an optically thick 
slab, of thickness L = b£, with 5 > 1 being large but finite. A generalization of the 
reasoning of section 2.1 allows us to write down the following Schwarzschild-Milne equation 
in this geometry 

r 6 (r, //, Ha) = po(A»> /Ua)e" r/Ma 

+ f dr' f%p Q {n^')e-^^'Y b {r'^'^ a ) 
Jo Jo l ^ 

+ J" dr' J 1 j^Po(^ -^')e- {T '-^^'Y h (r\ V, (2 _ 3g) 

+ fdr' / ,1 ^i2( / i / )po(/*,/* , )e- (T+T ' )/M 'r 6 (r , ,-/i / ^a) 
Jo Jo A ^ 

fdr' f ^R{-n')po{n, V)e- (26 - T - r ' )/M 'r6(TV,^). 
Jo Jo 
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The solution of this equation for a thick slab can be constructed from both solutions 
Ts and Th of the half-space geometry, by means of a matching procedure, along the lines 
of ref. [12]. Using the asymptotic forms (2.25), we obtain 

F s (t, ii, Ha) - rj^^ r H (r, n) (r finite, b - r » 1), 
r 6 (r,^,^ a )^<; ° + ZToT (2.40) 

& + 2 * Fg(& " T ' ~^ {b " T finit6 ' T 1} - 
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Both expressions lead to a linear (diffusive) behavior [3, 12] in the bulk of the sample 
(r»l, 6-r>l), namely 

r b(r, /*, lia) = ^~t~~ [b - r + ror* + ^(r* - 1)] + 0(e~\ e"^) . (2.41) 

The derivation of the diffuse transmitted intensity per solid angle element dO& again 
closely follows the lines of ref. [12]. We obtain 

dT( ?-" b) = T * A T (9 a ,9 b ), (2.42) 
dQ b b + 2r T* y ' ; ' v 1 

with 

A T (9 a ,9 b ) = ™ Sda J aTb TMTM, (2.43) 
127rm z \x a \x h 

where we have again used some of the definitions of Table 1: m = no/ni is the index 
mismatch, and T a = T(// a ) and T b = T(fj, b ) are the transmission coefficients in the incident 
and outgoing directions, respectively. 

The most non-trivial part of the above result consists in the function Ti(/x), whose 
definition and general properties have been exposed in section 2.2. It will be evaluated in 
the large index mismatch regime in section 3, and in the very anisotropic regime in section 
4. 

The result (2.42) shows that the effective thickness of the sample is (b + 2tot*)£ = 
L + 2tq£*- In other words, z$ = t$£* represents the thickness of a skin layer. This quantity 
is also referred to as the injection depth, or the extrapolation length of the problem. Our 
result thus shows on a firm basis that the relevant length scale of the problem in the 
diffusive regime is the transport mean free path £*, rather than the scattering mean free 
path £. Let us anticipate that such a simple picture is not always correct. Our expressions 
of the width of the backscattering cone and of the absorption length will indeed involve 
both mean free paths through the combination \J 11* . 

2.5. Enhanced backscattering cone 

The general formalism exposed above can be extended to the study of the celebrated 
enhanced backscattering phenomenon, which takes place in a narrow cone in the vicinity 
of the exact backscattering direction [4]. This phenomenon is due to the constructive 
interference between any path in the medium and its time- reversed counterpart. It is well- 
known that, for isotropic scattering, the width of the cone is of order A9 ~ Xi/£ ^ 1. One 
of the goals of this section is to derive a quantitative estimate of the width of the cone, 
with emphasis on its dependence on the anisotropy of the scattering mechanism. 

As recalled in section 1, the form of the backscattering cone can be predicted by 
resumming the so-called maximally-crossed diagrams. We restrict the analysis to normal 
incidence (9 a = 0), and to the geometry of a half-space diffusive medium. We introduce a 
dimensionless transverse wavevector 

Q = q£ = k £9' = h£9, (2.44) 
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where 9' and 9 are the incidence angles of the outgoing radiation, and ko and k\ are its 
wavenumbers, respectively inside and outside the diffusive medium. 

According to ref. [12], the reflected intensity in the vicinity of the backscattering 
direction, i.e., for 9 <C 1, k\i 3> 1, and Q fixed, takes the form 

A C (Q) « ^ [ 7 (1, 1) + lc{Q) - Po(l, -l)/2] . (2.45) 

In this expression, 7(1, 1) represents the background reflected intensity in the normal 
direction, in agreement with eq. (2.38). The second term 7c (Q) is the contribution of 
the interference between the sequences of any number (n > 1) of scattering events and 
their time-reversed counterparts. The subtracted third term is the contribution of the 
single-scattering events (n = 1), which are invariant under time inversion, and should not 
be counted twice. We define the enhancement factor 

B{Q) = A^O) = 1 + WJ) ' ( } 

It will turn out that the peak value of the interference contribution coincides with the 
background term, i.e., 7c(0) = 7(1, 1), so that the enhancement factor at the top of the 
cone, namely 

B(0) = 2 - Po(1, ~ 1) , (2.47) 
{) 2 7 (1,1)' 1 } 

nearly equals two, up to the small contribution of single-scattering events. 

From a more quantitative viewpoint [9, 12], we have 



/>oo 

1C {Q)= / dre- T r c (Q,T,-l), 
Jo 



(2.48) 



where the function Tc(Q, t, y) is a solution of a Q-dependent Schwarzschild-Milne equation 
of the form 

r C (<9, r, fi) = po(fji, l)e" T + (M * r c )(Q, r, ^), (2.49) 
with a Q-dependent kernel M = Mb + Ml, in analogy with eq. (2.23), namely 

M B (Q, r, /*, r', = (?0*'(r - r '))£2^ e -(r-TW Jf) _ r /| , 



Ml(Q, r, /*, r', = ^( V)^7p^^( V)e- (T+T,)/I/I Jo (Q(t + r')^F^l 



(2.50) 

where J is the Bessel function. 

The shape of the top of the backscattering cone, described by the small-Q behavior 
of ic{Q)-> is expected to be universal. It is indeed due to the contribution of long paths 
in the diffusive medium, along which the radiation undergoes many scattering events. On 
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the contrary, the tails of the cone, corresponding to a large reduced wavevector Q, only 
involve short sequences of two, three, etc. scattering events, and are therefore expected 
to depend on the details of the scattering mechanism. This is already apparent in the 
subtracted term in eqs. (2.46, 47), which involves the single-scattering cross-section in the 
direction of exact backscattering. We shall come back to this point in section 5. 

The universal small-Q behavior of the backscattering cone can be determined analyti- 
cally as follows. Deep in the bulk of the medium, i.e., for r 3> 1, the integral Schwarzschild- 
Milne equation (2.49) can be approximated by a differential equation, obtained by expand- 
ing the function Tc(Q, t, y) in powers of (r' — r). Keeping only the first two derivatives, 
we obtain the following diffusion-like equation 



F c (Q,t,^ = V 



T C (Q, r, fj) - /j^-r c (Q, r, n) 



d 2 o 2 1 (2 ' 51) 

V-n>rc(Q, T tf i)-¥-(i- y 2 )r c (Q, r, + • • 



The asymptotic behavior of the source function Fc(Q,t, //) for r ^> 1 and Q <C 1 can be 
derived by looking for a solution to eq. (2.51) in the form of the following expansion in 
Legendre polynomials 

T C (Q, r, fi) = e~ s ° T c tPM- (2-52) 

£>0 

Along the lines of Appendix A, and making use of eqs. (2.14) and (A. 4), we obtain the 
following recursion relations for the coefficients q, involving the coefficients wt of the 
expansion (2.13) of the phase function 

ct f, Q 2 \ ( ^ , l + i 

1 7T ) c e + s 0\ Wn 7 C ^-1 + nn I o C W 



vj t \ 2 J \2£- 1 2£ + 3 

f 2 f £(£ - 1) 2£ 2 + 2£-l , (£ + !)(£ + 2) 

+ { S ° + 2 J \(2£ - 1)(2£ - 3)° e ~ 2 + (2£ - 1)(2£ + 3)° e + (2£ + 3)(2£ + 5)° e+2 

(2.53) 

It is clear from the structure of these relations that, when Q is small, sq scales as 
<5, and the coefficients q ~ decay rapidly. Keeping this hierarchy in mind, and using 
zuq = 1, we expand the relations (2.53) for £ = and £ = 1 as follows 

= ^r-^co + + • • • , (— - l) Cl = s c + ■■■ (2.54) 

Using eq. (2.16), we obtain 

so^^L, ^«( T *_i)J^L (Q«l). (2.55) 



The next step also follows the lines of ref. [12]. The small-Q expansion of the Q- 
dependent kernels (2.50) only contains even powers of Q. As a consequence, the term 
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proportional to \Q\ in Tc(Q, t, y) is proportional to the homogeneous solution r#(r, y) of 
the ordinary (Q = 0) Schwarzschild-Milne equation. Putting everything together, we are 
left with the following estimates of the source function Fc(Q, t, /x). For Q <C 1 and fixed 
t we have 

r c (Q, r, /*) = r s (r, ^, 1) - ri (l)^Lr H (r, + 0(Q 2 ), (2.56) 



whereas for Q <C 1 and r>l simultaneously we get 



T C (Q, r, n) = Tl (l)e-^/V^ ^ + M (^(r* - 1) - t t*) + 0(Q 2 )) . (2.57) 

The universal peak of the backscattering cone is then evaluated by inserting the esti- 
mate (2.56) into eq. (2.48), using eqs. (2.29, 31). We thus obtain the following expression 

7c (q) = 7 (i,i)(i-M + 0(q*)), (258) 



where the width of the cone reads 



i.e., in physical units 



AQ = W 7^' ( 59) 



A6= 3 -^^=, (2.60) 

n(i) 2 fciv^F 



with k\ being the wavenumber of radiation outside the diffusive medium. The 1/ V ££* 
scaling law of the width of the backscattering cone is one of the most striking results of 
this paper. We shall come back to this point in detail in section 5. 

2.6. Extinction and absorption lengths 

Up to now we have assumed that the diffusive medium is conservative. This means 
that the light only experiences elastic collisions; there is neither absorption, nor inelastic 
scattering, implying the normalization (2.2) of the phase function. We now want to discuss 
briefly the case of a weakly absorbing diffusive medium, characterized by a non-trivial 
albedo a such that 1 - a « 1. In this case the diffuse intensity is expected to die off 
exponentially inside the medium, with a characteristic absorption length L a b s . 

The known expression [2, 3, 9] of the absorption length can easily be recovered by 
means of the formalism exposed in the previous section, in the case of general anisotropic 
scattering and in the regime of weak absorption. We shall actually determine the extinc- 
tion length of the more general problem defined by the Q-dependent Schwarzschild-Milne 
equation (2.49). We look for a slowly varying source function of the form (2.52), along the 
lines of section 2.5. The main difference is that we have now w = a ^ 1. We expand the 
recursion relations (2.53) for £ = and £ = 1 as follows 

/ ' -l) c = $2 °~ Q2 c + S -^ Cl + (—-l) Cl = SQCo + ... (2.61) 



zuq J 3 3 



16 



from which we obtain 

4 « (2 . 62) 

The Q-dependent extinction length reads 

it x 1/2 



WQ) ~ - ~ n2 ■ \ (Q « 1, 1 - a « 1). (2.63) 



The absorption length is obtained by setting Q = in the above expression, namely 

(l-a«l), (2.64) 



7 abs 



\ 1/2 



3(1 -a) 



in agreement with refs. [2, 3, 9]. 

3. LARGE INDEX MISMATCH REGIME 

In this section, we extend to the general case of anisotropic scattering the approach of 
ref. [12], which predicts the behavior of physical quantities in the regime where the optical 
indices uq and n\ of the diffusive medium and of the surroundings are very different from 
each other, i.e., when their ratio m = n^/ni is either very small, or very large. As already 
pointed out in ref. [12], important simplifications occur in these regimes of a large index 
mismatch, where the mean transmission coefficient of the boundaries of the medium is 
very small. To be more specific, radiation cannot enter the medium (respectively, leave 
the medium) in the limit m <C 1 (respectively, m >> 1), except at normal incidence. Ref. 
[9] already contains part of the novel results of this section. 

3.1. Diffuse reflection and transmission 

Along the lines of ref. [12], we evaluate the reflected and transmitted intensities in the 
large index mismatch regime by means of the following singular perturbative expansion of 
the Green's function Gs(t, fi, t', fi'), defined by eq. (2.26). 

The starting point consists in noticing the identity 

M l (t, ii, t', h') = R(-fi')M B (r, fi, -r', -//) (3-1) 

between both kernels defined in eq. (2.23). Using R(fi) = 1 — T{pt), we can recast eq. 
(2.26) as 

G s (t, t\ h') = po(ji, h')6(t - t') 

+ f°dr" f ^f[M B (r-r"^,0^")+M B (r + r"^,0,-^")} 

J ° J ~ 1 f3 2) 

xGs{t",h",t>,h>) y ' J 

- J™ dr" j X ^-T{n")M B {r, t", -lj")G s (t", , r' , p'). 
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In the limit of an infinitely strong index mismatch, i.e., for m = or m = oo, the 
transmission coefficient T(p) vanishes identically, so that the last integral of eq. (3.2), 
involving T(/i), is absent. It can be checked that the remaining terms only determine the 
Green's function up to an additive constant. This constant is only fixed by the integral 
involving T((i) in eq. (3.2). It can therefore be expected to diverge asm^O and m — > oo. 

In order to demonstrate this explicitly, we expand the Green's function according to 

G s (t, (jl, r', n') = C S + G (r, //, r', //) + •••, (3.3) 

with the hypothesis that Cs diverges, whereas G (t, /j,, t', n') remains finite, and the dots 
stand for terms which go to zero, as m — > or m — > oo. The finite part Go(t, /j,, t', //) 
obeys the following equation 



G (r, n, t', fjf) = po(n, /i')S(t - t') 

- 1 d/i" 



+ J™ dr" J ^L[M b (t- r", /*, 0, fi") + M B (r + r", /*, 0, V')] 



xG (t'V,tV) 
- 1 d/x", 



(3.4) 



C S jf° dr" f ^ ^-T{h")M b {t, ^ r", V), 



together with the consistency condition 
, f 1 d/i f 1 dfi' 



J dT J dr' J ^T^)M B (r + r'^,0,-^)Go(r',^,r'',^)=0, 



(3.5) 



derived along the lines of ref. [12]. 

The constant Cs of the expansion (3.3) can be derived by integrating eq. (3.4) over 
the variables < r < oo and — 1 < ji < 1. This yields 

Cs = (3.6) 
where T is the mean flux transmission coefficient 

T = * m =\ 0.7) 

I 3m 2 (m+ l) 2 V ~ ; 

This quantity assumes its maximum T = 1 in the absence of any index mismatch, i.e., for 
m = 1, and it vanishes in both cases of a large index mismatch, according to 

— (m«l), 

T ~ { o ( 3 - 8 ) 
(m > 1). 
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The asymptotic behavior in the limits m C 1 and m ^> 1 of the quantities pertaining 
to reflection and transmission is immediately obtained by replacing in eqs. (2.29-31) the 
Green's function Gs(t, /U, t', /jf) by its leading constant term C5. We thus obtain 

r « — , ri(/i) « — , 7(// ,W>) ~ r ■ (3.9) 

These predictions are identical to those derived in ref. [12], in the case of isotropic scat- 
tering. We have thus shown that the quantities which determine the diffuse reflected and 
transmitted intensities do not depend at all on the anisotropy of the cross-section in the 
large index mismatch limit. 

3.2. Enhanced backscattering cone 

The shape 7c (Q) of the cone of enhanced backscattering for a normal incidence can 
also be evaluated analytically in the regimes m < 1 or m > 1 of a large index mismatch. 
To do so, we recast the Q-dependent Schwarzschild-Milne equation (2.49) in a form similar 
to eq. (3.2), namely 

T c (Q,t, //) =p (^,l)e~ T 

+ f°dr' / ^f[M B (Q,T-T',»,0,lJ>') + M B (Q,T + T',»,0,-»')} 

Jo J - X 2 (3.10) 

xr c (Q,rV) 

- J™ dr' jf * ^7V)M B (Q, r, ^ r', -n')T c {Q, r', yi). 

By inserting the estimates (3.9) into the general result (2.59), we anticipate that the 
width of the cone is small in the large index mismatch regime, since it is proportional 
to the mean transmission T . Our present goal is to determine the scaling law of the 
backscattering cone jc(Q) m the scaling regime where both Q and T are simultaneously 
small. Under these circumstances, we look for a solution to eq. (3.10) of the form 

r c (Q,T^)^ lc (Q)e-W T/V ^, (3.11) 

where the extinction length of the exponential decay has been taken from eq. (2.57). 
We insert this form into eq. (3.10), and we integrate again both sides over the variables 
< t < 00 and — 1 < \i < 1. The integrals independent of the transmission T(p) can be 
performed explicitly, by means of the formula 

/_'^(i + a (i-/. a )r v> = ! ^»i-0 a /s (Q«D, (3.12) 

whereas the Q-dependence of the integral involving T{ji) can be neglected in the scaling 
regime. After some algebra we are left with the following result 

lc{Q) " + r/4 W< 1 - r < 1 )- < 3 ' 13 > 
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The small-Q expansion of this prediction reads 



7c«)4(l-^ + -) (3.14) 
The width of the cone therefore scales as 

3T 

AQ w TTt' ( 3 - 15 ^ 
4Vt* 

in agreement with the general formula (2.59), together with the results (3.9). 

4. VERY ANISOTROPIC SCATTERING 

In this section we investigate the regime where the scattering cross-section is very 
anisotropic, i.e., strongly peaked in a narrow cone of width rms <C 1 around the forward 
direction, with 2 ms = (O 2 ). We thus have 1 — vo\ !=s 2 ms /2 < 1, so that the anisotropy 
parameter r* ~ 2/0 2 ms is very large. The wavelength Ao of radiation in the medium, 
the scattering mean free path £, and the transport mean free path I* can therefore be 
considered as three independent length scales (Ao <C i <C £*), besides other characteristic 
lengths, such as the sample thickness L, and possibly the absorption length L abs . 

The interest in this very anisotropic regime is twofold. First, we shall show that the 
Schwarzschild-Milne problem in the absence of internal reflections is exactly solvable in 
this regime, just as it is for isotropic scattering, whereas the intermediate situation of a 
general anisotropy can only be treated numerically. Second, the dependence of physical 
quantities on anisotropy, with respect to the well-known isotropic case, can be expected to 
yield the largest effects in the very anisotropic regime. We shall come back to this point 
in section 5. 



4.1. The example of large spheres 

We first show that very anisotropic scattering can be realized experimentally. We 
consider the multiple scattering of light by large dielectric spheres, with radius a much 
larger than the wavelength Ao in the medium, i.e., with scale parameter k^a ^> 1, so that 
geometrical optics can be used. Furthermore we assume that the optical index ns of the 
spheres is very close to the mean index uq of the medium, i.e., 

— =m s = l + S s (ks|«l). (4-1) 
n 

It is expected that scattering mostly takes place in the forward direction in this regime. 

The study of the scattering cross-section of spheres is an old classical subject. The 
full solution was first derived by Mie in 1908. Ref. [25] provides an extensive overview of 
this field. The regime k$a ^> 1 and \Ss\ <C 1 still has to be split into several subcases, 
according to the value of the combination \5s\koa. This can be understood as follows. In 
the framework of geometrical optics we can distinguish between the diffracted light, which 
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is outgoing within an angle ©diffr ~ l/(^o a )? independent of 5s, and the refracted light, 
which is outgoing within an angle O re fr ~ \$s\, independent of k$a. 

From now on we concentrate our attention on the regime \5s\ <C 1, k$a ^> 1, and 
\5s\koa 3> 1. In this regime we have Odiffr <C © r efr <C 1. The cross-sections associated 
with each of the above processes asymptotically reads Odiffr ~ o"refr ~ ™ 2 , so that the 
total elastic cross-section a m 2na 2 is twice the geometrical one. This is the celebrated 
extinction paradox. We make the following approximations. We neglect the diffraction 
phenomenon by setting Odiffr = 0. We treat the refracted light according to the laws of 
geometrical optics, as illustrated in Figure 1. We neglect all the rays which are reflected 
at least once at the surface of the sphere, so that we only have to consider the refracted 
ray drawn on the Figure. The corresponding phase function reads 



4 sin/? cos (3 



dp 



de 



(4.2) 



sin 6 

with the notations of Figure 1, which also imply 

« -25 s tan (3. (4.3) 
Eqs. (4.2, 3) lead to the Lorentzian-squared scaling form [25] 

1 (\S 2 

»-< e > x w+k?- (44) 

The cross-section is thus strongly peaked in the forward direction, as anticipated. 

We now turn to the evaluation of the coefficient tt7i )re fr corresponding to refracted 
light. The integral (0 2 ) re fr = 2 £>refr(0) ©d6, with the phase function of eq. (4.4), 
is logarithmically divergent. A more careful treatment is therefore needed, which consists 
in using directly eq. (4.2), choosing u = sin 2 (3 as integration variable. We thus obtain 

1 4 



G7l,refr = 7^2 + ^2 \ ^du yj (1 - u) (m| - u) , (4.5a) 



with w max being the smaller of both numbers 1 and m|. The integral in eq. (4.5a) can be 
performed explicitly for both ms < 1 and ms > 1. In both cases we obtain 

w lreb ss 1 -25fln-^-, withA = 2e" 3/2 (\5 S \ < 1). (4.5b) 

The anisotropy parameter can be derived from the above estimate, not forgetting that 
Odiffr ~ °"refr ~ o/2 implies 1 — W\ = (1 — zui )Te f r )/2. To leading order for \5s\ ^ 1, this 
yields 

£* 1 



= J ^ A- ^ 
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It is therefore evident that large spheres with a weak dielectric contrast provide a phys- 
ical instance of very anisotropic scattering, as described in the beginning of this section, 
with the refraction angle O re f r ~ \8s\ ^ 1 playing the role of © rms , up to a logarithmic 
correction. We shall come back to this last point in section 5. 

4.2. Scaling limit of the Schwarzschild-Milne equation 

We now show that the general formalism of section 2 undergoes important simplifica- 
tions in the regime of very anisotropic scattering. These will allow us to derive analytical 
results for physical quantities, including the shape of the backscattering cone, in the ab- 
sence of index mismatch at the interface. This solution, to be described in sections 4.3 
and 4.4, therefore shows for the first time that both limiting cases of isotropic scattering 
and of very anisotropic scattering are nearly on the same footing as far as the existence of 
exact results is concerned. 

The simplifications in the regime of very anisotropic scattering can be understood in 
physical terms as follows. Consider the random sequence of scattering events experienced 
by a light ray. At every scattering event, the direction n of the ray is only modified by a 
slight amount of order © rms , with 9 r 2 ms 2/t*. The extremity of the vector n therefore 
performs a Brownian motion on the unit sphere, with a small angular diffusivity e ~ ©rms- 
A similar picture has been used for long in the context of semiflexible polymer chains: the 
so-called Kratky-Porod theory models polymers as continuous persistent walks, whose unit 
tangent vectors obey a diffusion equation on the sphere (see ref. [26] for a review). 

In more quantitative terms, the relationship (2.7) between the n-dependence of the 
specific intensity I(r, n) and the source function T(r, n) takes the form 

r(r,n)w(l + eA)/(r,n), (4.7) 

where A is the Legendre operator (Laplace operator on the unit sphere) 

d d 2 Id 2 

A - cot0 a» + W + ^roW (4 ' 8) 

As a consequence, the linear operator V introduced in eq. (2.11) reads 

V^l + sA , (4.9) 

where 

is the Legendre operator (4.8) in the ^-independent sector. 

The angular diffusivity e is then determined by observing that P\ (fx) = \i is an eigen- 
function of the operator A , with eigenvalue (—2), and of the operator V, with eigenvalue 
zui, by virtue of eq. (2.14). We thus identify vo\ = \ — 2e, whence, using eq. (2.16), 

1 t 2 
e = — = — w^, (4.11) 

2r* 2£* 4 ' v ' 
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in agreement with the above heuristic estimate. 

The radiative transfer equation (2.8) thus assumes the differential form 



2t*h—I(t, fi) = A I(T, fx), 



(4.12) 



which demonstrates that physical quantities vary over length scales of order r ~ r*, i.e., z ~ 
£*. In the next two sections we present our exact analytical predictions concerning physical 
observables pertaining to optically thick slabs, in the regime of very anisotropic scattering, 
in the absence of internal reflections. The analysis roughly follows the presentation of 
section III of ref. [12], devoted to the case of isotropic scattering. 

4.3. Exact treatment in the absence of internal reflections: homogeneous 
Schwarzschild-Milne equation 

This section is devoted to an analytic determination of the solution r#(r, ix) of the 
homogeneous Schwarzschild-Milne equation in the regime of very anisotropic scattering, 
and of the associated quantities of interest, To and t\{lx). 

It is convenient to consider the Laplace transform gjj(s,Lx) of r#(r, ix), defined as 
follows 



drr H (r, -At)e ST (Res<0). 



9h(s,lx) = 

The Schwarzschild-Milne equation (2.24) can be recast as 

~9h(s,lx) 



(4.13) 



9h(s,h) 




1 + six 
9h(s,h) - T*n(ii)/3 



(-K/x<0), 
(0 <//<!), 



(4.14) 



1 + six 

since we have gu{—\j fx, fx) = t*t\(ix)/3 for < fx < 1, as a consequence of eqs. (2.31, 35). 

We now expand eq. (4.14), using the expression (4.9, 11) of the operator V. We 
introduce the rescaled Laplace variable 



a = 2sr*, 

as well as a new unknown function hu (cr, fx) , defined as follows 



9H{s,fx) 



4r* 2 (l + sfx)h H (a,fx) 



(-Kfj,<0), 



4t* 2 (1 + sfx)h H (o, fx) + t*ti(ai)/3 (0 < fx < 1). 



(4.15) 



(4.16) 



We assume that Ti(ix) vanishes faster than linearly in fx for fx — > 0. This hypothesis will 
be checked a posteriori, as it will turn out that Ti(fx) ~ fx 3 / 2 for fx — > (see eq. (4.53)). 
The function Iih(&, fx) is then continuously differentiable as a function of fx, and it obeys 
the following equation 



(A - aix)h H (cr,ix) 



(-1<a*<0), 
ti(a0/6 (0 <//<!). 



(4.17) 
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By means of the change of variable 

fx = tanhx, (4-18) 

eq. (4.17) can be recast in the form of the following inhomogeneous Schrodinger equation 
on the real x-line 

hr B (.,x)-<,v(x)h H (.,x) = {° v{xHx)/e [*<g; (4.i9) 

In this formula the accents denote differentiations with respect to x. The potential 

V(x) = tanhx(l - tanh 2 x) (4.20) 
is an odd function of x, such that V(x)dx = fidfx, and we have set 

T 1 (ji)=fiu(x) (0 < At < 1, x>0). (4.21) 

Equation (4.19) is a self-consistent equation for the two unknown functions hn(o;x) 
and v(x). Analyticity properties in the complex a- variable turn out to allow for an exact 
analytical solution of this equation. 

We introduce a basis of two elementary solutions {u±(a, x), 112(0; x)} of the (homoge- 
neous) Schrodinger equation 

u"(a,x) - aV(x)u(a,x) = 0, (4.22) 

with the following asymptotic behavior 

ui(a, x) ~ 1, U2(cr, x) ~ x (x — > — 00), (4.23) 

up to exponentially small corrections. Similarly we introduce the basis {vi(a, x), V2(cr, x)}, 
with the asymptotic behavior 

vi(a,x)pal, V2((T,x)ttx (x — > 00). (4.24) 

The Schrodinger equation (4.22) admits solutions with the above boundary conditions, 
since the potential V(x) vanishes exponentially as x — > ±00. The four above solutions are 
entire functions of a, i.e., they are analytic in the whole cx-plane. Furthermore they are 
related by the following identities 

v 1 (a,x) = Ui(-a, -x), v 2 (a, x) = -u 2 (-cr, —x), (4-25) 

because the potential V(x) is odd. 

We recall that, if u(x) and v(x) are any two solutions of the Schrodinger equation 
(4.22), with the same a, their Wronskian 

W{u, v} = u(x)v'(x) - u'(x)v(x) (4.26) 
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is independent of x. The boundary conditions (4.23, 24) imply that both bases of functions 
have unit Wronskian, namely 

W{ Ul ,u 2 } = W{v u v 2 } = 1. (4.27) 

For generic values of the parameter a, both bases of solutions are related by a 2 x 2 
transfer matrix of the form 

vi)-{h(*) F(-a)){u 2 )- (4 ' 28) 

The three functions which enter eq. (4.28) are entire functions of a; the determinant of 
the matrix is F(a)F(—a) —G(a)H(a) = 1; as a consequence of eq. (4.25), G(a) and H(a) 
are even functions of a. 

The above functions also govern the non-trivial asymptotic behavior of the bases of 
solutions 



u\(a, x) m F(— a) — G(a)x, u 2 (a, x) m —H(a) + F(a)x (x — > oo), 
vi(a, x) a* F(a) + G(a)x, v 2 ((T, x) « H(a) + F(— a)x (x — > — oo), 

as well as their mixed Wronskians 



(4.29) 



W{u 1 ,v 1 } = G(a), W{ Ul ,v 2 } = F(-a), W {u 2 ,v 1 } = -F(a), W{u 2 ,v 2 } = -H(a). 

(4.30) 

The first terms of the Taylor expansion of ui(a, x) and u 2 (o~, x) around a = read 



Ui(a, x) = 1 — — (1 + tanhx) + • • • , 

u 2 (cr, x) = x + a ^— (1 — tanhx) + ln(2 cosh a:) j + 



(4.31) 



We have also determined the terms of order cr 2 , which are too lengthy expressions to be 
reported here. They imply 

F(a) = 1 + a + a 2 /2 + ■■■, G(a) = a 2 /3 + ■■■, H(a) = (7/6 - n 2 /36)a 2 + ■■■ (4.32) 

On the other hand, large values of the complex parameter a correspond to the semi- 
classical regime for the Schrodinger equation (4.22), where the behavior of the functions 
u±(a,x) and u 2 (a, x) can be derived by means of a W.K.B.-like approximation. This 
regime is analyzed in detail in Appendix B. Let us mention that the wavefunctions display 
oscillations when either a > and x < or a < and x > 0, whereas they are growing or 
decaying exponentially in the other two cases. 

The function G(a) deserves some more attention, since it will play a central role in the 
following. G(cr) can be viewed as the functional determinant of the Schrodinger equation 
(4.22), in the following sense. Assume a is such that G(a) = 0. We have then 

v\{<j,x) = F{a)ux{a,x), m(a, x) = F{-a)v x {a, x), F(a)F(-a) = 1. (4.33) 
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In other words, for such a a, the Schrodinger equation (4.22) has a bounded solution 
over the whole real line. Throughout the following, using slightly improper terms, such 
a value of a is called an eigenvalue of the Schrodinger equation, and the corresponding 
function vi(a,x) is referred to as the associated eigenf unction. The semi-classical analysis 
of Appendix B demonstrates that there is an infinite sequence of real eigenvalues. We label 
them by an integer — oo < n < oo, so that a n > for n > 1, o~o = 0, and cr_ n = — o~ n . The 
associated eigenfunctions v\(a n ,x) are orthogonal with respect to the following indefinite 
metric 

/•OO 

vi(a m ,x)vi(a n ,x)V(x)dx = N n 5 m ^ n (-oo < m, n < oo). (4.34) 



/ 

J — i 



The squared norms N n obey the symmetry property 



N- n N n 



F(-o- n ) F(a n ) 



(4.35) 



as a consequence of eq. (4.33). The eigenfunction i>i(0, x) = 1 associated with the zero 
mode o"o = is peculiar, since its squared norm reads Nq = 0. As a consequence, the basis 
of eigenfunctions {vi(a n ,x),n ^ 0} only spans the set of bounded functions f(x) on the 
real x-line such that 

/oo 
f(x)V(x)dx = 0. (4.36) 
-oo 

For such functions, we have 

/oo 
v 1 (a n ,x)f(x)V(x)dx. (4.37) 

The determination of the contribution of the zero mode to functions f(x) which do not 
obey the condition (4.36), such as e.g. a constant, requires more care. An elegant way of 
dealing with this problem consists in introducing a deformation parameter as we shall 
see in section 4.4. 

It is advantageous to factor the entire function G(a) as follows 

G(a) = ^P(a)P(-a), (4.38) 



with the notation 



p(v) = n f 1 + f ) ■ ( 4 - 39 ) 



The explicit solution of the inhomogeneous Schrodinger equation (4.19) now goes as 
follows. Since hu{o-,x) is a regular function of x in the x — > — oo limit, it is proportional 
to Wi(cr, x) for x < 0, namely 

hn{o, x) = an{(y)ui{a, x) (x < 0). (4.40) 
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For x > we solve eq. (4.19) by varying the constants, namely we look for a solution of 
the form 

Jih(o~, x) = &h(ct, x)vi(a, x) + ch(o~, x)v2(a, x) (x > 0). (4-41) 
The unknown constants bn(cr, x) and ch(ot, x) obey the requirements 

b' H (a, x)v 1 (a, x) + c' H (a, x)v 2 (a, x) = 0, (4.42a) 
b' H (a,x)v[(a,x) +c' H (a,x)v' 2 (a,x) = V(x)v(x)/6, (4.42b) 

where the accents again denote differentiations with respect to x. Eq. (4.42a) is a con- 
straint imposed a priori, in order to break the large redundance of the representation 
(4.41); eq. (4.42b) is then a consequence of eq. (4.19). Eq. (4.42) can be solved for b' H 
and c' H , using eq. (4.27). Integrating the expressions thus obtained, we get 



poo 

b H (a, x) = b H (o; oo) + / v 2 ((J, y)u(y)V(y)dy/Q, 

J X 

/>oo 

c H (o-,x) = - v 1 (a,y)v(y)V(y)dy/6. 

J X 



(4.43) 



Indeed there cannot be a non-zero constant ch(o~, oo), because this would correspond to 
an unacceptable singular solution of the form hn{o, x) ~ x ~ ln(l — y) for fj, — > 1. 

Both expressions (4.40) and (4.41) have to match at x = 0, together with their first 
derivatives. These two conditions determine 6h(ct, 0) and c#(c, 0), and some algebra then 
leads to the identity 

/•OO 

G{a)a H {a) = -c H (a,0) = / v 1 (a,x)u(x)V(x)dx/6. (4.44) 

Jo 

We can argue on eq. (4.44) as follows. Since gH(s,fi) is analytic in the half-plane 
Res < 0, au((y) has the same property for Re a < 0, so that the zeroes — a n of G{a) 
cannot be poles of afl(a). They are therefore zeroes of c#(cr, 0). On the other hand, the 
a n 's are not zeroes of ch{o~, 0), since the integral expression in eq. (4.44) is positive for 
a > 0. Hence they are poles of a#(o"). The final step concerns the small-cr behavior of the 
quantities involved in eq. (4.44). The asymptotic behavior (2.25) of Th(t,h) yields the 
following double-pole structure for its Laplace transform near s = 

g H (s,ri = \- T * {T0 + ^ +0(1) ( a ->0), (4.45) 
s z s 

which implies the following small-cr behavior of an(o") 



We thus obtain finally 



a H {a) = -1 + 1-22 + 0(1) (<t-0). (4.46) 
la 



a ^) = a 2 P \_ a y -ch(*,0) = ^. (4.47) 
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Eq. (4.47) can be considered as an explicit result. Indeed -P(cr), denned in eq. (4.39), is for 
all purposes a known function, since the eigenvalues a n can be determined numerically, es- 
sentially with arbitrary accuracy, via the partial- wave expansion procedure of Appendix A. 
Furthermore the semi-classical analysis of Appendix B determines the asymptotic behavior 
of the eigenvalues in the regime of large quantum numbers (n >> 1). 

As a first consequence of the above results, we can determine the reduced thickness To 
of a skin layer, i.e., the reduced extrapolation length, by comparing the small-cr behavior 
of the expression (4.47) for a#(a) with the expansion (4.46). We thus obtain 

r = 1 - 2 V — = 0.71821164. (4.48) 

n>l n 

The series converges, since a n grows as n 2 , according to the semi-classical estimate (B.21). 
The number given in eq. (4.48), as well as all the subsequent ones, has been obtained by 
means of the partial- wave expansion described in Appendix A. This number gives an idea 
of the accuracy of this approach. The most significant numerical results are listed in Table 
2, together with their counterparts in the case of isotropic scattering. 

Second, the determination of ti(/j,) goes as follows. We recall that this quantity is 
related to u(x) by eq. (4.21). Eqs. (4.44, 47) yield 

/•OO 

/ v 1 (a,x)u(x)V(x)dx = 2P(a). (4.49) 
Jo 



The small-a expansion (4.31) of u\(a, x), together with eq. (4.25), allows us to recover 
the sum rules (2.34, 37b) in the following form 

/•OO /-OO 

/ v{x)V{x)dx = 2, / u(x)tanhxV(x)dx = 2t . (4.50) 
Jo Jo 



On the other hand, the semi-classical estimates (B.18, 20) of vi(a,x) and P(cr) for 
large values of a = K 2 imply 

/•oo 

/ u(x)Ai{K 2 ' 3 x)xdx « (G/n^K- 5 / 3 (K > 1). (4.51) 
Jo 

This estimate yields, by means of eq. (B.25) for s = 3/2, the small- a; behavior of u(x), i.e., 

u{x) w 6(2x/tt) 1/2 (x«1). (4.52) 



We thus obtain the following universal scaling behavior of the Ti-function in the case 
of very anisotropic scattering 

r 1 (^)^6(2/7r) 1 /V /2 (//<!)• (4-53) 
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This novel result is in contrast with the linear behavior Ti(/i) ~ ji\/?> observed for isotropic 
scattering. The law (4.53) confirms the hypothesis made in the beginning of this section, 
namely that Ti(p) vanishes faster than linearly in \i. 

Finally, we can extract the full functions v{x) and T\(\i) from eq. (4.49) by means of 
the inversion formula (4.37). We obtain 

v{x) = 3(r + tanhx) + 2 ^T^iK, x), (4.54) 

n>i iVn 

i.e., 

^ = 3(r + iu) + 2 4^^iK, argtanh^). (4.55) 

n>l 

The semi-classical analysis of Appendix B implies that the contribution of the n- 
th mode to the results (4.54, 55) falls off exponentially with n for x > 0, according to 
exp [ — K n {l — 7(x))]. This exponential convergence disappears as x — > 0, where eqs. 
(4.52, 53) apply. The explicit terms in front of the sums in eqs. (4.54, 55), corresponding 
to the contribution of the zero mode (n = 0), have been anticipated from results to be 
derived in section 4.4. 

Figure 2 shows a plot of the function T±(p), for both isotropic scattering [12] and very 
anisotropic scattering (eq. (4.55)). The maximal values ti(1) for both cases are given in 
Table 2. The difference between both limiting cases is remarkably small. 

4.4. Exact treatment in the absence of internal reflections: inhomogeneous 
Schwarzschild-Milne equation and backscattering cone 

In this section we derive analytical expressions in the regime of very anisotropic 
scattering of the special solution Ts(r, n, fi a ), with its by-product the bistatic coefficient 
tG^o) /•*&); an d of the source function Tc(Q, t, fx), yielding the shape 7c (Q) of the enhanced 
backscattering cone at normal incidence. 

It is advantageous to view the above quantities as two different particular cases of 
the following more general problem, which has no direct physical interpretation. We con- 
sider the Q-dependent source function Fc(Q, t, /j,, /J, a ), which obeys the inhomogeneous 
Q-dependent Schwarzschild-Milne equation (2.49), with the source term of eq. (2.22), and 
we define the associated Q-dependent bistatic coefficient as 



poo 

lc(Q,Ha,Hb) = I dTe~ T/tlb rc(Q,T,-n b ,n a 

JO 



). (4.56) 



Notations are consistent, in the sense that we have 7( / u a , / u&) = 7c(<3 = 0, // a ,//b) and 
lc{Q) = IciQiHa = A*6 = 1). Among other advantages, this approach will treat in an 
elegant way the problem of the zero mode with vanishing norm of the Schrodinger equation 
(4.22). 
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In analogy with eq. (4.13), we define the Laplace transform of the Q-dependent source 
function as follows 



poo 

gc(Q,s,fi,/j,a) = / drr c ((5, r, -fi,fi a )e ST (Res < 0). 
Jo 



(4.57) 



The Schwarzschild-Milne equation (2.49) can be recast as 



g c (Q,s,ji,ii a ) 



with the definition 



^aPo(-^,Mo) 
W(Q, S, -fl a ) 



+ v 



9c(Q,s,n,n 
w(Q,s,n) 



V 



w(Q,s,n) 



(-K/x<0), 
(0 <//<!), 



(4.58) 
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(4.59) 



In analogy with section 4.3, we expand eq. (4.58), using eqs. (4.9, 11). We introduce 
the rescaled variables 

(t = 2st\ k = Q\^* = hVife, (4.60) 



assuming k > 0, as well as a new unknown function hc(K, a, /jl, /j, a ), defined as follows 

(-K/x<0), 



(n , _ j 2T*w(Q,s,n)h c (K,(r,fi,fia) 



+ 7c(Q,A»,A»a) (0 <//<!). 



(4.61) 



The definition (4.60) of the rescaled variable k is in full accord with the general 
prediction (2.59) for the width of the backscattering cone. In more precise terms, and in 
analogy with eq. (2.58), we have 

7c(Q,Ha,Hb) = n/(p a ,l*b) - gTi(// a )ri(/i b ) + 0(k 2 ) (k«1). (4.62) 

We assume that 7c(Q 5 A* a , /■*&) vanishes faster than linearly as \i a or ^6 — > 0. Using the 
change of variable (4.18), we are again left with an inhomogeneous Schrodinger equation, 
namely 



h f Q(K, cr, x, x a ) — (aV(x)+K 2 W(x))hc(K, cr, x, x a ) 



-2/i a 8(x + x a ) (x<0), 



(4.63) 



fj, a V(x)p(K,x,x a ) (x>0). 
The potential has, besides the odd component V(x) of eq. (4.20), an even component 

W(x) = (1 -tanh 2 x) 2 , (4.64) 

and we have set 



lc{Qi l^i H>a) = WaP(K, x a ) (/i = t&nh x > , n a = tanhx a > 0). 



(4.65) 
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The main difference introduced by the deformation parameter k resides in the spec- 
trum of the Schrodinger equation 

u"(a, x) — (aV(x) + k 2 W{x))u{o,x) = 0. (4.66) 

The eigenvalues with label n^O acquire a regular ^-dependence of the form 

<r n («) = -<t_„(k) =a n + 0(n 2 ), (4.67) 

as well as the associated eigenfunctions v±(k, cr n (K), x), whereas the double degeneracy of 
the zero mode o"o = is lifted into the following two exact eigenvalues and eigenfunctions 
ofeq. (4.66) 

<7 ±0 (k) = ±2k, ui(k, (t ±0 (k), x) = e ±K ( 1 - tanhx ) = e ±K(1 "' i) , (4.68) 
with squared norms 

W ±0 (/c) = ±^ " cosh 2k) . (4.69) 

The introduction of labels ±0 is consistent with setting n = in formulas such as eq. 
(4.35) or (4.67), rather than with the standard arithmetics of integers! 

For any k ^ 0, the set of eigenfunctions {vi(k, c n (/c), x)}, where the label n runs over 
the non-zero algebraic integers (n ^ 0) plus both values n = ±0, now spans the whole 
space of bounded functions f(x) on the real line. The difficulty of the constraint (4.36), 
due to the vanishing norm of the zero mode at k = 0, is thus cured in a natural way. 

The mixed Wronskian G(n,a) = W{ui,vi} can still be factorized over its zeroes, in 
analogy with eq. (4.38), namely 

G(k, a) = G(k, 0)R(k, <t)R(k, -a), (4.70) 

with 



n>0 



The prefactor G(k, 0) of eq. (4.70) is a non-trivial function of k. Indeed this quantity 
can be viewed as the functional determinant of the Schrodinger equation (4.66) with a = 0, 
namely 

u"(cr,x) - K 2 W(x)u(a,x) = 0. (4.72) 

Eq. (4.72) is equivalent to the spheroidal equation studied by Meixner and Schafke [27]. 
Some properties of this equation have been studied in detail [28], in an investigation of 
nematic phases of semiflexible polymer chains. The occurrence of eq. (4.72) in that context 
is related to the Kratky-Porod description of persistent chains, mentioned in section 4.2. 
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Eq. (4.72) has a discrete spectrum of imaginary eigenvalues of the form k — i*£ n 
(n > 1), as shown by the semi-classical analysis of Appendix B. On the other hand, the 
regularity of G(k, a) at k = implies the small- k behavior G(k, 0) ~ An 2 /3, whence 

G^0) = ^l[(l + ^). (4.73) 

6 n>l V *„/ 

The solution of the inhomogeneous Schrodinger equation (4.63) follows the lines of 
section 4.3. The particular values x = —x a < and x = define three sectors, in which 
we look for a solution of the following form 

{a c (K, a, x a )ui(n, a, x) (x < -x a ), 

b c (K, a, £ a )wi(K, a, x) + c c (k, a, x a )u 2 (K, a, x) (-x a < x < 0), 

d c (K, a, x, x a )vi(n, a, x) + e c («, a, x, x a )v 2 (K, a, x) (x > 0). 

(4.74) 

The constants which enter the last of these expressions obey the conditions 



d' c (K, a, x, Xajv^K, a, x) + e' c (K, a, x, x a )v 2 (K, cr, x) = 0, 

d! c (K, a, x, x a )v[(K, a, x) + e' c (K, a, x, x a )v 2 (K, a, x) = p a V(x)p(K 7 x, x a ), 



(4.75) 



whence 



poo 

dc{K, cr, x,x a ) = dc{K, a, oo,x a ) + p a / v 2 (k, cr, y)p(K, y, x a )V(y)dy, 

J x 

POO 

e c( K ? °~i x i x a) = -p a \ v 1 {K,a,y)p{K,y,x a )V{y)dy. 

J x 

On the other hand, the matching of the solution (4.74) at x = —x a yields 

&<?(«> cr, x a) = ac(«i 0", ^o) + 2VclU 2 (k, a, x a ), 
c c {k, a, x a ) = -2/j, a ui(K, cr, x a ), 

whereas its matching at x = leads to 

cr, 0, x a ) = F(k, -a)b c (K, a, x a ) - H(k, ct)c c (k, cr, x a ), 
-e c («, cr, 0, x a ) = cr)6 c («, cr, x a ) - F(«, o-)c c (k, cr, x a ), 



(4.76) 



(4.77) 



(4.78) 



with the notations (4.28). We are thus left with 

G(k, (j)ac(K, cr, X a ) = -2/IaV^K, (7, -x a ) - € C (k, (7, 0, x a ). (4.79) 

We can now follow the approach used on eq. (4.44). Since cic(k, a, x a ) is holomorphic 
in the half-plane Re a < 0, the zeroes —a n {K) for n > of G{n,a) cannot be poles of 
ac(«, cr, x a ). Hence they are zeroes of the right-hand side of eq. (4.79). 
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We are therefore left with the problem of finding 6c(k, a, 0, x a ), an entire function of 
a, from the knowledge of its values on the sequence of points a = —a n (K) (n > 0), together 
with a natural assumption of minimal growth at infinity compatible with these data. This 
is a generalization to an entire function of the problem of finding a polynomial Q(z) with 
minimal degree, knowing its values at N points, namely Q(z n ) = Q n for 1 < n < N. It is 
useful to view the z n 's as the zeroes of the normalized polynomial 

P(z)= H (z-z n ). (4.80) 

l<n<7V 

The solution Q(z) with minimal degree (generically N — 1) is given by the following La- 
grange interpolation formula 

«w= E «. n r=r- = PW E r^TWirr (481) 

l<n<N l<m^tn<N l<n<N 

Extending eq. (4.81) to the present case of an infinite sequence of data for an unknown 
entire function, we obtain 

-e c (K, a, 0, x a ) = 2» a R( K , a) £ ^ wTp/h^ V ^ < W ^ 

f-^ Q (a + a n {K))(dR/da){K, -<r n {K)) 

or equivalently, using a generalization of the identity (C.6) to k ^ 0, together with the 
definition (4.70), 

-e c («, o-, 0, x a ) = -2/j, a G{K, 0)R{k, a) > — ■ r^^TT~\ • ( 4 - 83 ) 



n>0 



Finally, we can derive an explicit expression for p(K,x,x a ), by means of an inversion 
formula analogous to eq. (4.37), namely 



I \ nnt n\ ^(«»^m(«))-R(«,^n(«)) Vi(k, <T m (K) , x) V 1 (k, a n (n) , X a ) 

p{K,X,X a ) = -2G(/C,0) > rT - r—r -r 

™7>0 ^m(«)+a n (K) N m {Kj N n {K) 

_ ^ Yl^ll ^n(K),x)vi(K, Q- n («Q, -X ) + l?i («, (7 n («) , X a )t>l («, CT n («),-x) 



iV n («) 

n>0 

(4.84) 

This central result is manifestly symmetric under the exchange of x and x a , as it should. 
This symmetry property is a valuable check of the whole approach, since both arguments 
x and x a have played uneven roles throughout the derivation. 

We now consider the n — > limit of the above results, which will then be denoted with 
an S'-subscript instead of a C-subscript, for the sake of consistency with the notations of 
section 2. Eq. (4.83) can be recast for k — > as 

( n ^ o p/ Ji a ^ a n P(a n )vi((T n ,x a ) \ 
-e s ((7, 0, x a ) = 2p a P{a) 1--^ ^— ttt • (4-85) 

V (^ + ^)iv n y 
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First, we are now able to complete the proof of the anticipated result (4.54, 55), by inserting 
eq. (4.85) into eq. (4.79), expanding the latter equation for k — > to the first non-trivial 
order as a — > 0, and comparing the result with the estimate 

a s (a,x a )*-^^- ((t-0). (4.86) 
a 

Second, the small-cr expansion of eq. (4.85) allows us to recover the sum rules (2.33, 37a) 
in the following form 

poo poo 

/ p{x, x a ) V(x)dx = 2, / p(x, x a ) tanhxy(x)da: = 2v(x a )/3 — 2 tanhav (4.87) 
Jo Jo 

On the other hand, for large values of a = K 2 , we can use the semi-classical estimates 
(B.18, 20, 22), setting a n = /c 2 , in order to transform eqs. (4.76, 85) into 

/•oo poo t.2/3 

I p(x,x a )Ai(K 2 / 3 x)xdx^(2/n)K 1/3 — 2 — M(k 2 / 3 x a )dk . (4.88) 

o Jo k + K 

This estimate shows that xx a p(x, x a ) is a homogeneous function of its arguments with 
degree zero, when both of them are small, p(x,x a ) = g(x/x a ). The rescaling of 

eq. (4.88) according to z = K 2 l 3 x = k 2 l 3 x a then yields, by means of a mere identification 
of both integrands, using eq. (B.25), the expression g(z) = (3/n)z 3 / 2 /(z 3 + 1), implying 
the following scaling behavior 

3(x a x b ) 1/2 

p(x a , x b ) « a — 3- (x a , x b < 1), (4.89) 



/ 



n{xl + x 3 ) 



or equivalent ly 



7(fc,^)^ _rf , GW6«1)- (4.90) 



3(^ b ) 3 / 2 

This novel result is in contrast with the rational behavior 7(/z a , /Life) ~ p a Pb/(Pa + Pb) 
in the case of isotropic scattering. The law (4.90) confirms the hypothesis made in the 
beginning of this section, namely that 7(/U a , p b ) vanishes faster than linearly in either of its 
arguments. It is also worth noticing that the scaling form (4.90) of the bistatic coefficient 
saturates the sum rule (2.33). 

The full expression of the bistatic coefficient 7(/U a , p b ) is obtained by taking the k — > 
limit of eq. (4.84), using the definition (4.65). The modes m, n = yield divergent 
contributions as k — > 0, which cancel out as they should, as well as finite parts, so that we 
are left with the result 

7 ^ a '^ = 3r + ^(p a + Pb) + 2 ^2 P j^ n [vi(<r n ,x a ) +vi(a n ,x b )] 

PaPb * n>1 

- XI 7v~ [ v i( a n, x a )vi(a n , -x b ) + vi(a n , -x a )vi(a n , x b )] (4.91) 

n>l n 

2 a m a n P(a m ) P(a n ) 
~q TV j^M ff m,Xa)Vl{<T n ,X b ), 

lib - lb ^ _1_ 
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with /i a = tanha: a , fit, = tanhx^. 

The maximal value of the bistatic coefficient, which yields an absolute prediction for 
the diffuse reflected intensity in the normal direction, reads 7(1, 1) = 4.889703. This 
number is some 15% above the corresponding one in the case of isotropic scattering (see 
Table 2). 

The shape 7c(«) of the enhanced backscattering cone at normal incidence in the 
regime of very anisotropic scattering can also be derived from the general result (4.84). 
Indeed, by letting \i a = [i^ = 1 in the latter formula, and using again a generalization of 
the identity (C.6), we obtain the following exact result 

M = 2V - 2C(k 0) V R(K,a m (K))R(K, g TC («)) 

1C[ } ^ (dG/d*)(«, a n (K)) [ ' } J^ Q (a m ( K ) + a n {K))N m {K)N n {n) • 

(4.92) 

For small values of the reduced wavevector, the following behavior 

7c(«)=7(l,l)-r 1 (l)V3 + C(« 2 ) («<1) (4.93) 

is more directly obtained from eq. (4.62). On the other hand, for large k, the second 
(double) sum in eq. (4.92) is negligible, whereas the first (simple) one can be evaluated 
by means of the semi-classical estimate (B.35). After some algebra using the identity 
lim z ^i J2n>o(~ z ) n = 1/2, we are left with the estimate 

7c(k) « 2ne~ 2K (k > 1). (4.94) 



Figure 3 shows a plot of the enhancement factor B(Q), defined in eq. (2.46), against 
the reduced wavevector k = Qy/r*, for both isotropic and very anisotropic scattering. 

A last comment is in order. For a large but finite anisotropy, the shape of the enhanced 
backscattering cone 7c (Q) only follows the universal scaling law (4.92) for k of order unity, 
i.e., Q of order 1/y/r* ~ 6 rms . For larger values of Q, the backscattering cone exhibits 
non-universal tails, which are expected to be generically exponentially small. 

4.5. Extinction lengths of azimuthal excitations 

Up to this point, we have only considered samples in the form of a slab, or a half- 
space, and investigated physical quantities which have the cylindrical symmetry around 
the normal to the slab, i.e., which are independent of the azimuthal angle <p. 

We now want to consider briefly the </>dependence of the diffuse intensity in the regime 
of very anisotropic scattering. Because the sample has cylindrical symmetry, the azimuthal 
dependence of physical quantities, such as the diffuse intensity 7(r, n) = 7(r, 9, ip), can be 
expanded over the modes e zrnip , where the algebraic integer m is the azimuthal, or magnetic 
quantum number. Different sectors associated with different values of m are decoupled. 
As already mentioned in section 2, all the sectors contribute e.g. to the reflected intensity, 
except if the incident beam is normal to the sample, or, more generally, has itself cylindrical 
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symmetry. The situation is different in transmission through thick slabs, to which only 
the sector m = 0, considered so far, contributes. The reason for this is that the intensity 
in the other sectors is exponentially damped inside the sample, namely 

J m ~e- 2 /S (4.95) 

where l m is the extinction length of the azimuthal excitations in the sector ra- 
it is the purpose of this section to determine these lengths in the limit of very 
anisotropic scattering. The Legendre operator (4.8) reads 

d <9 2 ra 

in the sector defined by the azimuthal integer ra. As a consequence, and along the lines of 
sections 4.3 and 4.4, we are led to study the Schrodinger equation 

u"(a, x) — (ra 2 + aV(x))u{a,x) = 0. (4.97) 

The corresponding extinction length is given by 

2£* 

Zm = — , (4.98) 

^0 

where is the smallest positive eigenvalue of eq. (4.97). This is indeed the location of 
the first singularity of the Laplace transform of the intensity. 

For large values of the azimuthal integer ra, we can use the semi-classical analysis 
developed in Appendix B. The Sommerfeld quantization formula associated with eq. (4.97) 
reads 

J {-aV(x)- ra 2 ) 1/2 dx = J i^Y^ ~ (1 ^ 2)2 j « (n + 1/2)tt. (4.99) 

The smallest eigenvalue corresponds to setting n = in the above formula. In first 
approximation we express that the argument of the square-root inside the //-integral in eq. 
(4.99) has zero as its maximal value, so that [i- = We thus obtain ~ 3y / 3ra 2 /2. 
In second approximation we expand the integrand around its maximum, which takes place 
for fi ~ y/3/3. We obtain after some algebra the following next-to- leading order estimate 

a ( m) w (3^/2) (ra 2 + mV2), (4.100) 

whence 

— = =r- (ra>l). (4.101) 

3v/3(ra 2 + my/2) ' ' 
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This semi-classical estimate gives accurate numbers of the whole spectrum of extinc- 
tion lengths, down to the largest one, t\. Indeed eq. (4.101) yields ~ 0.318861, 
whereas the exact numerical value reads = 0.2829169. 

For a large but finite anisotropy, the lengths l m follow the universal law (4.101) only 
for m < m* ~ \fr* ~ 1/G rms . For larger values of the azimuthal number, the i m become 
non-universal numbers of order I. This crossover is expected on physical grounds. Indeed 
large azimuthal numbers m ^> m* correspond to an angular resolution 50 <C O rms , so that 
the details of the cross-section matter in this regime. 

5. DISCUSSION 

In this paper we have considered several aspects of the theory of multiple scattering 
of scalar waves, with a general anisotropic cross-section. We have limited this study to the 
geometry of an optically thick slab, of thickness L » I*. Two different length scales are 
relevant, namely the scattering mean free path I and the transport mean free path £*. The 
general results of section 2 concern the dependence of various quantities on both mean free 
paths. Quantities for which only the long-distance diffusive behavior matters scale with 
the transport mean free path £*. This is the case for e.g. the angle- resolved transmission 
through a thick slab (2.42), and especially for the thickness z = tqI* of a skin layer, i.e., 
the extrapolation length. It will indeed be recalled later on that the anisotropy of the 
scattering mechanism only contributes a small effect, entirely contained in prefactors such 
as the constant To, or the functions Ti(fx) and 7( / u, / u / ). 

We have also obtained an unexpected result concerning the width of the backscat- 
tering cone, scaling as A9 ~ Xo/VIi* (2.60). This prediction is in disagreement with the 
previously accepted scaling law in Xq/£* [29, 30]. The derivations of the latter result rely 
on the diffusion approximation. A reason why this approximation may fail will be given in 
a while. It is clearly apparent from sections 2.5, 6 that the derivation of our scaling law for 
the width of the backscattering cone is fully similar to that of the known law (2.64) of the 
absorption length [2, 3], namely L a bs ~ Vli*. This connection between the backscattering 
peak and the absorption length is seemingly rather accidental. As a matter of fact, such 
a relationship is to some extent already present in the diffusion approximation [9] . Indeed 
in this framework all extinction effects, provided they are small enough, can be coded in a 
single parameter, namely the mass M. such that 

M 2 = q 2 + zfi + (5.1) 

abs 

where q = Q/£ is the transverse wavevector, L abs is the absorption length, and O = 
(oj — a/)/-D p hys represents the properly dimensioned contribution of a small frequency shift 
between the advanced and the retarded amplitude propagators which build up the diffuson. 
The scaling laws (2.60) and (2.64) for the width of the cone and the absorption length can 
be put together as M 2 ~ 1/ {it). 

In order to understand in heuristic terms the role played by the intermediate length 
scale in the case of the absorption length, we can use the following argument, based 

on the random- walk picture of diffusive propagation. If the absorption cross-section is very 
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small, i.e., if the albedo is very close to unity, a photon will only be absorbed after a large 
number of collisions, of order N ~ 1/(1 — a). The absorption length is the typical depth 
reached by a photon during N steps of its random walk through the diffusive medium. 
In the regime of very anisotropic scattering, the random walk is persistent over a length 
t > £. We thus have L^ bs ~ N*£* 2 , where N* is the effective number of independent 
steps of length t. The total path length reads Nt ~ N*£*, whence N* ~ £/[£*{! - a)], 
and L^bs ~ ££*/(! — a). This is precisely eq. (2.64), up to a numerical factor. 

The role of the intermediate length scale y/Ii* in the shape of the backscattering cone 
is more difficult to understand at an intuitive level. A photon entering the diffusive medium 
at the origin will leave the medium at a random transverse position . In rough terms, the 
width of the cone scales as A9 ~ Ao/Ar^, where Ar± is the typical transversal extent |rj_| 
of the path of the photon in the medium. In the case of isotropic scattering, the well-known 
result AQ ~ 1 corresponds to Ar± ~ £. The paths which determine the width of the cone 
thus have a transversal extent of order one mean free path. As a consequence, even though 
the triangular shape of the top of the cone is determined by the longest paths, which are 
Brownian, and thus correctly described by the diffusion approximation [29, 30], the width 
of the cone is determined by the shortest paths, and thus not necessarily quantitatively 
described by the diffusion approximation. If we extend the above line of reasoning to 
the regime of very anisotropic scattering, we are led to consider paths with the shortest 
transversal extent, such as that shown symbolically on Figure 4. Each step makes an angle 
comparable to the maximal allowed angle O rms with respect to the previous one, so that 
the path consists of N ~ 1/ rms steps. Its total length and its transversal extent scale as 
Ar± ~ N£ ~ v '11* . This explains the scaling law for the width of the cone, provided such 
stressed paths have the correct length scale. The statistics of the transversal extent r± of 
persistent random walks has been the subject of recent works [31, 32]. Unfortunately the 
results of numerical simulations presented there only concern the universal diffusive fall-off 
of the distribution for r± 3> £* , so that the above heuristic idea can neither be confirmed 
nor infirmed. 

The next question of interest concerns to what extent physical quantities are universal, 
and to what extent they depend on the particular microscopic scattering mechanism. As 
far as the diffuse reflected or transmitted intensities are concerned, the detailed structure 
of the scattering mechanism only contributes a small effect, contained in prefactors of the 
laws mentioned above, such as the constant To, or the functions Ti(//) and 7( / u, / u / ). Two 
situations allow for a more detailed investigation of these matters. 

(i) The regime of a large index mismatch, where the boundaries of the sample almost act 
as perfect mirrors, is considered in section 3. Our main result (3.9) shows that the physical 
quantities considered do not depend at all on the scattering cross-section in this regime. 
This can be easily understood as follows. Since the thickness z$ ~ 4£*/(3T) of a skin layer 
is very large, the radiation undergoes many scattering events near the boundaries before 
it leaves the medium, so that the details of every single scattering event are washed out. 

(ii) In the absence of internal reflections, we have considered in detail the regime of very 
anisotropic scattering. Our exact treatment, based on the expansion (4.7), is expected to 
be valid in the regime O rms «C 1 of a broad universality class of phase functions. Although 
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this universality class cannot be easily characterized, we can assert that it contains at least 
the phase functions scaling as 

P (e)^$(e/e rms ), (5.2) 

such that the scaling function $ has a finite second moment. This restrictive definition 
does not encompass a priori the Lorentzian-squared phase function (4.4), which has a 
logarithmically divergent second moment, as already mentioned in section 4.1. The same 
remark holds for the so-called Henyey-Greenstein phase function 

P{G) = (1-2, cose + ^)3/^ (5 - 3) 

often used in numerical investigations [3, 14], for which the second-moment integral is lin- 
early divergent. In section 4 we have presented an exact analytical treatment of universal 
aspects of very anisotropic scattering, yielding to the results (4.48, 55, 91, 92) for the 
observables pertaining to thick slabs. These results are not entirely explicit, since they are 
given in terms of the eigenvalues and eigenfunctions of one-dimensional Schrodinger equa- 
tions, which are accessible both numerically, via the partial-wave expansion of Appendix 
A, and analytically in the limit of large quantum numbers, via the semi-classical analysis 
of Appendix B. The exact solution derived there is nevertheless an analogue in the regime 
of very anisotropic scattering of the exact solution in the case of isotropic scattering, which 
dates back to Chandrasekhar [1] . 

The discussion of the dependence of physical quantities on the details of the scattering 
mechanism is summarized in Table 2, where we compare the numerical values of the dimen- 
sionless absolute prefactors of five characteristic quantities, for isotropic scattering and for 
very anisotropic scattering. The relative differences, shown in the last row, are very small 
in most cases. Some other quantities, such as the shape of the enhanced backscattering 
cone, or the spectrum of extinction lengths of the azimuthal excitations, exhibit universal 
behavior in the very anisotropic regime only in a limited range, corresponding to a low 
enough angular resolution (50 3> © rm s), so that the details of the scattering cross-section 
do not matter. 

Finally, we can compare our universal results in the very anisotropic scattering regime, 
for some of the quantities listed in Table 2, with the outcomes of numerical approaches. Van 
de Hulst [3, 14] has investigated in a systematic way the dependence of various quantities 
on anisotropy, for several commonly used phenomenological forms of the phase function, 
including especially the Henyey-Greenstein phase function (5.3). The data on the skin- 
layer thickness reported in ref. [3] show that, as a function of anisotropy, To varies from 
0.7104 (isotropic scattering) to 0.7150 (moderate anisotropy), passing a minimum of 0.7092 
(weak anisotropy). The trend shown by these data suggests that our universal value 
0.718211 is actually an absolute upper bound for To- Numerical data concerning ti(1) 
is also available. Van de Hulst [14] has extrapolated two series of data, concerning the 
Henyey-Greenstein phase function (5.3), which admit a common limit for very anisotropic 
scattering (g — > 1). According to the analysis of section 4, this limit reads in our language 
n(l)/4 = 1.284645, whereas ref. [14] gives the two slightly different estimates 1.273±0.002 
and 1.274 ± 0.007. The agreement is satisfactory, although it cannot be entirely excluded 
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that the observed 0.8% relative difference can be a small but genuine nonuniversality effect. 
Indeed, as mentioned above, the Henyey-Greenstein phase function (5.3) might not belong 
to the universality class where our approach holds true. The same remark applies to a 
less complete set of data [14] concerning the intensity 7(1, 1) of reflected light at normal 
incidence. 
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Appendix A. Partial-wave expansions 

In this Appendix we describe a numerical algorithm based on a partial- wave expansion, 
that we have used to determine the eigenvalues and eigenfunctions of the Schrodinger 
equations involved in section 4. 

We first consider the Schrodinger equation (4.22). Going back to the variable, this 
equation reads 

(A -<7//)vi(<7,//) = 0, (A.l) 

where Ao is the Legendre operator in the ^-independent sector, defined in eq. (4.10). It 
is natural to expand the function v\ (a, fx) in the Legendre polynomials 

v 1 (<T,n) = ^2a e ((r)P e (fj). (A.2) 

£>0 

Indeed these polynomials are eigenfunctions of Ao, namely 

A P e = -£(£+l)P e , (A.3) 
and the product fJ>Pe(n) has the following expression 

(2£ + l)fiPM = (£+l)P e+1 (iu)+£P e -M, (A.4) 
so that eq. (A.l) amounts to the following three-term recursion relation 

£(£ + l)a t + a (J±L ai+1 + -^—^^ = 0. (A.5) 

When a is one of the eigenvalues a n , eq. (A.5) has an acceptable solution {ag(a)}, 
decaying to zero for large £. The quantities needed in section 4 can then be evaluated as 
follows. The normalization condition (4.24) becomes 



^2a t (<7 n ) = l, (A.6) 



since Pe(l) = 1. The squared norms N n of the eigenfunctions read 

N » = 4 E (2 e+\)(L + 3) ae{an)ae+lM (A - 7) 

as a consequence of the normalization of the Legendre polynomials 



^PkMPM = (A.8) 
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Finally, the non-trivial mixed Wronskians F(±a n ) read 

F(-a n ) = = Vl (a n , fji=-l) = £(-1) e a e (c n ), (A.9) 

since P e (-1) = 

We now consider the ^-dependent Schrodinger equation (4.66), namely 

(A -<Tfj, + k 2 (/U 2 - 1))ui(k, cr, At) = 0. (A. 10) 

We again expand the wavefunction over the Legendre polynomials 

ui(k,ct,//) = ^2a £ (K,a)P e (iJ,). (A.ll) 

By iterating eq. (A. 4) twice, we obtain the following five-term recursion 



- K " {(2i- l)(2i-3) ae - 2 + (2i-l)(2i + 3) ae + WTWl+t) ae+2 ) ~ °- 

(A.12) 

Eqs. (A.6, 7, 9) still hold true. 

We finally consider the wave equation 

(A-<TfJL)v 1 (<T,fJL,<p) = 0, (A.13) 

where A is the full Legendre operator, defined in eq. (4.8). Since the potential does not 
involve the azimuthal angle (p explicitly, we look for a solution v\ proportional to e im< ^, 
with m > being an integer. It is now natural to expand the function vi(a, fM,(p) in the 
spherical harmonic functions, which we choose to write as follows 

vifo /*,¥>) = e im * a e ,m(<T)Pr(»), (A-14) 

£>m 

in terms of the Legendre functions P™(//), which obey 

A(Pp(/j)e imLp ) = -£{£+l)Pp(fj)e imip , (A.15) 
and the product /jPp(fi) has the following expression 

(21 + l)fiPP(fi) = (£ + 1- m)Pp + M + (£ + m)Pp_M, (A.16) 
so that eq. (A.13) amounts to the three-term recursion relation 

(£ _|_ m _|_ \ £ _ m \ 

The recursion equations (A. 5, 12, 17) are easily implemented numerically. 



£{£-!) 2(1 -£-£ 2 ) t (£+!)(£ + 2) 
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Appendix B. Semi-classical analysis 

The outcomes of the semi-classical analysis presented in this Appendix are used at 
various places in section 4. We first consider the Schrodinger equation (4.22). For large 
values of the complex parameter cr, we look for rapidly varying solutions of the form 

u « $(x)~ 1/2 expj $(y)dy, (B.l) 

with 

$( x ) 2 = aV(x). (B.2) 

This approach is analogous to the celebrated W.K.B. approximation in Quantum Mechan- 
ics, since the condition a ^> 1 is equivalent to % being small. 

Because the potential V(x) is odd, we can restrict the analysis to the domain Re a > 0. 
We introduce the notation 

K = y/a. (B.3) 
Eq. (B.2) has real solutions for x > 0. We set 

p(x) = y/V(x) (x > 0), (B.4) 

so that $>(x) = Kp(x). Eq. (B.l) thus yields the basis of functions 

1 

(Kp(x)Y 
with 

/■OO 

I{x) = / p(y)dy (x > 0). (B.6) 

J X 

The above functions u± (x) are exponentially blowing up or decaying. The domains x > 
and a > (and similarly x < and a < 0) are said to be classically forbidden. 

On the other hand, for x < 0, eq. (B.2) has imaginary solutions. We set 

q(x) = V-V(x) (x<0), (B.7) 
so that <&(x) = iKq(x). Eq. (B.l) thus yields the basis of functions 

u± (x)w ] -— 7 =exp(±iKI(x)), (B.8) 



u±{x) w — / ^ 1/2 exp ( ± KI(x)), (B.5) 



with 



I(x)= j q(y)dy (x < 0). (B.9) 

«/ — OO 
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The above functions u±{x) are oscillating and bounded, up to the prefactor in q(x)~ 1 ^ 2 . 
The domains x < and a > (and similarly x > and a < 0) are said to be classically 
allowed. 

The difficulty of the semi-classical analysis comes from the existence of three turning 
points, namely x = and x — > ±oo, where the momentum variable p(x) or q(x) vanishes. 
The estimates (B.5, 8) loose their meaning in the vicinity of the turning points, where a 
more careful analysis is required, to be presented now. 

We first investigate the basis of functions {^1,112} for x < 0. For x — > —00 and 
K — > 00, the Schrodinger equation (4.22) assumes the simpler form 



u" + (2Ke x ) 2 u = 0. 



(B.10) 



A basis of solutions to this equation is given by the Bessel functions Jq(z) and Nq(z), with 
z = 2Ke x being a scaling variable. The boundary conditions (4.23) have to match the 
known small- z behavior of the Bessel functions, whence 



ui (a;) « J (2Ke x ), 

u 2 {x) « (n/2)N (2Ke x ) -QnK + lE )J {2Ke x ) (x 



-00), 



(B.11) 



where 7^ denotes Euler's constant. The known large-z behavior of the Bessel functions 
fixes the amplitudes of the integrals in eq. (B.8) for m and 112, namely 



u x (x) 



ivKq(x) 



1/2 



cos (KI(x) - 7r/4), 



/ 2 \ 1/2 

u 2 (x) « ( ^^J [(V 2 ) sin (if/(x) - tt/4) - (In if + 7b ) cos (K/(x) - tt/4)] . 

(B.12) 

On the other hand, for a; — * and K" — > 00, the Schrodinger equation (4.22) assumes the 
simpler form 



u" + i^:ru = 0, 



(B.13) 



which is equivalent to Airy's equation. A basis of solutions is given by the Airy functions 
Ai(z) and Bi(z), with z = K 2 l 3 x being again a scaling variable. The known behavior for 
z — > —00 of the Airy functions has to match eq. (B.12), whence 

u^x) « 2 1 / 2 K~ 1 / 3 [sm(KI) Ai(K 2 / 3 x) + cos(KI) Bi(K 2 / 3 x) 
u 2 {x) « 2 1/2 K" 1/3 {(tt/2) - cos(KJ) Ai(K 2/3 x) + sin(KI) Bi(K 2/3 x) (B.14) 



-(In if + 1e ) \sm(KI) Ai(K 2 / 3 x) + cos(KI) Bi(K 2 / 3 



x 



for x — > 0, with 
7 = 1(0) = 



p(x)da; 



d// 
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1/2 



¥ T(3/4) 
2 r(5/4) 



1.198140, 



(B.15) 



this definition being consistent with eqs. (B.6, 9). 

We now investigate in a similar way the functions {vi, V2} for x > 0. For x — > 00 and 
K — > 00, a basis of solutions is given by the modified Bessel functions iio(z) an d ^o(^)) with 
z = 2Ke~ x . The boundary conditions (4.24) have to match the known small- z behavior 
of the Bessel functions, whence 

Vl ps I (2Ke~ x ), 

V ; B.16 

v 2 ~K {2Ke- x ) + {\nK + lE )I {2Ke- x ) (x -> 00). 

The known large- z behavior of the Bessel functions only fixes the amplitude of the solution 
u+ of eq. (B.5) for v\ and V2, namely 



/ 1 \V 2 

^ ~ V27rtfp(a0j eK/(X) ' (B.17) 
^(z) ~ (InK + 7b)ui(x) (x > 0). 



Finally, the known behavior of the Airy functions as z — > 00 yields 

^ 1 (a:)^2 1 /2K-i/3 Ai ( K 2/3 x)e KJ 5 
v 2 (x) (InK + 7 B )vi(x) (x->0). 



(B.18) 



The above expressions (B.14, 18) of both bases of solutions as x — > and if — > 00 
allow us to derive the following semi-classical estimates for the elements of the transfer 
matrix introduced in eq. (4.28) 

F(a) ps [sm(KI) - (2/n)(lnK + 1e ) cos(KI)]e KI , 
G(a) ~ -(2/tt) cos(K/)e KJ , 
if(cr) ps (lnK + 7j5 )F((j), 
F(-a) ps (lnK + 7j5 )G((T) (a = if 2 ^00). 

The estimate for G(cr) directly yields the following expressions for its factors P(±cr), 
defined in eq. (4.39) 



P{a) ps (3/tt) 



ki 

!/ 2 _ _ 



K " ' (B.20) 
P(-a) ps -2(3A)^E^) (a = if 2 ^00). 

We also obtain from eq. (B.19) an estimate of the eigenvalues a n = if 2 , which are 
the zeroes of G(a), in the form 

K n ps (n + 1/2) j (n>l). (B.21) 
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This semi- classical formula gives a very accurate description of the whole spectrum of the 
Schrodinger equation (4.22). Indeed the relative error is maximal for the first non-zero 
eigenvalue, for which eq. (B.21) predicts K\ ps 3.933086, i.e., some 3.2% above the exact 
numerical value K\ = 3.811562. 

The semi-classical expression of the squared norms N n can be evaluated by inserting 
the estimates (B.19) into the identity (C.6). We thus obtain 

N n ps 1 e ^ n i ( w > !)_ ( B 22 ) 

nK n 

In section 4 we also need the expression of the Mellin transform of the Airy function 
Ai(x), namely 

/•OO 

m(s)= / x s Ai(x)dx (Res>-1), (B.23) 
Jo 

which we have not found in standard handbooks. The Airy equation implies the functional 
equation 

m(s + 3) = (s + l)(s + 2)m(s), (B.24) 
whose correctly normalized solution is 

m « = 3 " ,/3 TW5)- (B ' 25) 

We now generalize the above approach to the ^-dependent Schrodinger equation (4.66). 
In a first step we only consider the spheroidal equation (4.72), corresponding to a = 0. The 
semi-classical regime corresponds to |«| ^> 1; the following two cases have to be considered 
separately. 

First, for k real, the analogues of eqs. (B.ll, 12) read 

ui(k, x) ~ Io(2Ke 2x ) (k — > oo, x — > — oo), 

Ul(K,x)w e -(l+tanhx) / < q n (B.26) 

(7r K (l-tanh 2 x)) 1/2 

Since the potential W(x) is even, we have v\(k,x) = u\{k, —x), whence 

2 

G(k,0) = 2ui(«,0)wi(K,0) « -e 2K (k -> oo). (B.27) 

7T 

Second, for k imaginary, we set k = i£. The analogues of eqs. (B.ll, 12) read 
ui(k, x) Jo(2^e 2x ) (£ — > oo, x — > — oo), 



- . ™o _L fon^ ^ _ *rlA\ ( a- ^ (\\ (B.28) 

(tt£(1 -tanh 2 x))" 



(^(1 + tanhx) -7r/4) (x < 0) 
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The spheroidal equation (4.72) has an eigenvalue whenever either ui(k, x) or its derivative 
with respect to x vanishes for x = 0. We thus obtain eigenvalues of the form 
well as the semi-classical estimate 

i n « (n + 1/2)tt/2 (n>l). (B.29) 

In the general case where both a and k are non-zero, only the spectrum of eq. (4.66) 
will be needed. Hence we can content ourselves with the Sommerfeld quantization formula. 
A similar treatment is used in section 4.5 for the full Legendre operator. 

For the Schrodinger equation (4.22) the Sommerfeld formula reads 

/ ( - aV(x)) 1/2 dx ^(n + 1/2) tt, (B.30) 

J — oo 

yielding again the estimate a ±K 2 , with if n as in eq. (B.21). 

For the ^-dependent Schrodinger equation (4.66), the Sommerfeld formula reads 

£ + ( - aV(x) - k 2 W{x)) 1/2 dx = (°Y^ ~ ^) ' dfi ~ (n + 1/2) 7r ' (B - 31) 

the integral being extended over the classically allowed domain, where the square-root is 
real. 

This implicit equation (B.31) for the semi-classical estimate of the eigenvalues ±a n (K) 
can be investigated in several limiting cases of interest. For small values of k, the integrand 
can be expanded in a straightforward way. We thus obtain 

aM 1 / 2 » K n (n) = (n + 1/2) ^1 + ^+1/2)* + " ') ^ ^ K <<C 1)- (B ' 32) 

This expression confirms the general result (4.67). On the other hand, for large k, the 
allowed domain shrinks down to a small vicinity of \i = 1. We obtain the spectrum of an 
effective harmonic oscillator, namely 

a n (n) w 2(2n+ 1)k (k > 1). (B.33) 

This formula holds for all values of n, down to n = 0. The definition (4.71) thus yields 



Finally, putting together the estimates (B.27) and (B.34), and using the definition (4.70), 
we obtain the following asymptotic formula 

G(k,<t) w -e 2K cos^ («»1). (B.35) 

7T 4« 



47 



Appendix C. Useful identities on the mixed Wronskians 

In this Appendix, we derive the identity (C.6) used in section 4, and more generally 
we give alternative expressions for the derivatives with respect to the spectral variable a 
of the mixed Wronskians F(a) : G{a), and H{a), which enter the transfer matrix (4.28). 

To do so, we start by considering the derivatives 

U a {a,x) = — (a = 1,2), (C.l) 

which obey the inhomogeneous Schrodinger equation 

U'J l (a,x)-aV(x)U a (a,x) = V(x)u ce (a,x) (a = 1,2). (C.2) 

These equations can be solved explicitly by varying the constants, along the lines of sections 
4.3 and 4.4. We thus obtain 

/X fX 
u 1 (a,y)u 2 (a,y)V(y)dy + u 2 {<j,x) J ul(a,y)V(y)dy, 

U 2 (a,x) = -ui(a,x) / u 2 2 {a,y)V{y)dy + u 2 {a,x) \ u 1 (a,y)u 2 (a,y)V(y)dy. 

J — oo J — oo 

By taking the x — > oo limit of the above expressions, and using the asymptotic be- 
havior (4.29), we get the following expressions 



dF(a) 

da 
dG(a) 

da 
dH{a) 

da 
dF(-a) 

da 

with the definition 



G(a)N 22 (a) + F(a)N 12 (a), 
-F(a)N 11 (a)-G(a)N 12 (a), 
-F(a)N 11 (a)-G(a)N 12 (a), 
-F(-a)N 12 (a) - H(a)N 11 (a), 



(C.4) 



/oo 

u ce (a,x)u[3(a,x)V(x)dx («,/?= 1,2). (C.5) 
-oo 

On the spectrum, i.e., for a = a n , we have G(a n ) = 0, by definition. Furthermore eq. 
(4.33) implies iVn(cr n ) = N n /F 2 (a n ), whence the identity 

dG(a)\ N n 

= -N n F(-a n ), (C.6) 



J a=an F(a n ) 

with N n being the squared norm of the eigenfunction vi(a n , x), defined in eq. (4.34). This 
result is very general; it also holds for the ^-dependent Schrodinger equation (4.66). 
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CAPTIONS FOR TABLES AND FIGURES 



Figure 1: Laws of geometrical optics for the refraction of light by a large dielectric sphere. 

Figure 2: Plot of exact expressions for T\{ji) in the absence of internal reflections. Dashed 
line: isotropic scattering, after ref. [12]. Full line: very anisotropic scattering (this work). 

Figure 3: Plot of exact expressions for the enhancement factor B(Q) of backscattering 
cone at normal incidence, in the absence of internal reflections, against k = Qy/r* = 
k\\flt*d. Dashed line: isotropic scattering, after ref. [12]. Full line: very anisotropic 
scattering (this work). 

Figure 4: Schematic plot of a stressed path yielding the shortest value of the transversal 
extent r± compatible with the opening angle © rms . 

Table 1: Conventions and notations for kinematic and other useful quantities. 

Table 2: Comparison of the numerical values of various quantities of interest from the 
exact solutions in the absence of internal reflections. First row: isotropic scattering, after 
ref. [12]; second row: very anisotropic scattering (this work). tq£* is the thickness of a 
skin layer; ti(1) and 7(1, 1) respectively yield the transmitted and reflected intensities in 
the normal direction; B(0) is the peak value of the enhancement factor at the top of the 
backscattering cone; Ak is the width of the backscattering cone, in units of k = Qy/r* = 
kiy/Ii*9. The third row gives the relative difference of the second case with respect to the 
first one. 
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Table 1 





outside medium 


inside medium 


optical index 


n\ 


no = mni 


wavenumber 


k\ = niu/c = 27r/Ai 


ko = uquj/c = mk\ = 2-7r/Ao 


incidence angle 


9 


9' 


parallel wavevector 


p = ki cos 9 
= ko a//U 2 — 1 + 1/m 2 


P = k cos 9' 
= 


total reflection 
condition 


m < 1 and sin 9 > m 
(i.e. P imaginary) 


m > 1 and sin6>' > 1/m 
(i.e. p imaginary) 



transverse 
wavevector 


|q| 


= q = k\ sin 9 = ko sin 9' = ko a/1 — « 2 


azimuthal angle 




reflection 
and 
transmission 
coefficients 


partial 
reflection 

total 
reflection 


R _(P-pV _ U-^H 2 -l + l/mA 
\P + p) \ f ,+ y^-l + l/m 2 J 

4Pp 4^a/^ 2 - 1 + 1/m 2 
(P 1 p) (/x + vV " 1 + I/™ 2 ) 

[ R=l 

T = 
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Table 2 





TO 


n(l) 


7(1,1) 


5(0) 


Ak 


isotropic 


0.710446 


5.036475 


4.227681 


1.881732 


1/2 


very anisotropic 


0.718211 


5.138580 


4.889703 


2 


0.555543 


A(%) 


1.1 


2.0 


15.7 


6.3 


11.1 
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